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CHEMEQ2:  A  Solver  for  the  Stiff  Ordinary  Differential 
Equations  of  Chemical  Kinetics 


1  Introduction 


This  report  documents  CHEMEQ2,  the  latest  version  of  the  FORTRAN  package  CHEMEQ  [1,2].  The 
CHEMEQ2  routines  integrate  sets  of  coupled,  nonlinear  ordinary  differential  equations  (ODEs)  of  the  form 


dyi 


1  <  2  <  n, 


(1) 


where  t/,*  is  the  density  of  the  species  and  gi  is  its  rate  of  change.  Our  primary  application  of  Eq.(l)  is  to 
sets  of  coupled,  nonlinear  ODEs  that  represent  chemical  reaction  sets.  In  this  case,  the  dependent  variables 
{y,}  are  concentrations  or  densities  of  reacting  chemical  species.  Sometimes  this  equation  is  supplemented  by 
another  equation  for  the  change  in  temperature  or  energy  release  that  results  from  the  species’  interactions. 
The  source  term  yi,  which  is  a  function  of  the  concentrations  and  the  thermodynamic  state,  may  be  written  as 
the  difference  of  the  production  rate  qi  and  the  loss  rate  Piyi.  The  timescales  Tj  =  1/pi  for  the  various  species 
differ  by  many  orders  of  magnitude  and  there  may  be  strong  coupling  between  species  (i.e.,  the  Jacobian 
matrix  dgi/dyj  has  significant  off-diagonal  elements).  Under  these  circumstances,  the  set  of  equations 
represented  by  Eq.  (1)  is  considered  stiff  and  does  not  lend  itself  readily  to  numerical  solution  by  classical 
methods  such  as  the  low-order  Euler  methods  or  higher-order  Adams-Moulton  methods  [3-5],  Such  a  system 
then  requires  special  techniques  to  obtain  an  accurate  solution  efficiently. 

The  coupled  reaction  set  represented  by  Eq.  (1)  is  often  a  part  of  a  larger  model  that  solves  these  equations 
coupled  to  the  partial  differential  equations  describing  fluid  dynamics.  In  such  cases,  chemical  reactions  are 
only  one  of  several  processes  that  might,  for  example,  include  advection,  diffusion,  or  radiation  transport. 
The  numerical  methods  commonly  used  to  solve  such  chemically  reacting  flows  use  process  splitting  (or 
operator  splitting)  [5].  The  basic  idea  in  operator  splitting  is  to  calculate  the  effects  of  individual  physical 
processes  separately  for  a  chosen  global  timestep  A/^,  and  then  combine  the  results  in  some  way.  Each 
process  in  turn  can  change  different  system  variables  during  Aig.  Then,  when  it  is  time  to  integrate  the 
ODEs  representing  the  chemical  changes  during  Atg,  the  integrator  is  presented  with  a  new  initial  value 
Manuscripl  approved  May  9,  2001 . 

I 


problem  in  each  computational  cell.  The  integrator  must  therefore  solve 

/.Ox  n 

~  9i>  yi{t  )  —  Pi  ^2) 

tot  =  t°  +  At,.  The  ODE  integration  may  subdivide  At,  into  smaller  steps,  At,  to  obtain  an  accurate, 
stable  solution.  Here,  the  timestep  At  is  called  the  chemical  timestep  because  it  is  the  timestep  that  the 

ODE  integrator  uses  to  advance  the  chemical  reactions.  The  size  of  At  generally  varies  during  the  course  of 
the  calculation. 

Given  that  fluid  dynamic  calculations  are  seldom  accurate  to  better  than  a  few  percent,  any  requirement 
of  the  chemical  integrator  to  calculate  the  species  concentrations  more  accurately  than  a  few  tenths  of  a 
percent  is  usually  excessive.  Therefore,  the  chemical  integrator  may  be  relatively  low-order.  Also,  since 
the  integrator  must  solve  multiple  initial  value  problems  “from  scratch”  at  every  global  timestep.  it  is 
necessary  to  use  a  single-point  method,  requiring  information  only  from  the  current  time  level  to  calculate 
the  concentrations  at  At.  This  is  in  contrast  to  multi-point  methods  that  must  store  concentration  or  source- 
term  values  from  several  successive  timesteps  in  order  to  advance  the  solution.  Within  the  calculation  for  a 
given  At„  multi-point  methods  have  a  start-up  penalty  until  a  sufficient  number  of  steps  have  been  taken  to 
build  the  history  required  for  the  calculation,  and  they  often  require  interpolation  procedures  if  At  changes 
during  the  integration.  The  chemistry  integration  from  the  previous  At,  does  not  provide  the  history  needed 
to  restart  the  integration  because  the  fluid  dynamics  calculation  changes  the  state.  The  values  at  the  end  of 
the  previous  chemistry  integration  are  therefore  not  the  values  at  the  start  of  the  next  chemistry  integration. 

By  comparison,  a  single-point  method  has  minimal  start-up  penalty  at  the  beginning  of  an  integration  step 
and  there  is  never  a  fluid  dynamic  inconsistency. 

CHEMEQ  [1,2]  is  a  second-order  single-step  ODE  integrator  that  has  been  used  successfully  as  a  part 
of  a  number  of  different  types  of  reacting-flow  codes.  These  have  included  applications  to  combustion  [6-10] 
and  solar  physics  [11-13].  CHEMEQ  is  a  hybnd  method,  which  means  it  chooses  between  a  stiff  method 
and  a  non-stiff  method  for  integrating  each  ODE  within  the  system  depending  upon  the  timescale  of  that 
equation.  CHEMEQ  has  been  shown  to  outperform  standard  stiff  ODE  solvers  by  a  factor  of  50-100  in  speed 
in  validation  studies  on  chemical  integrations  alone  (i.e.,  not  coupled  to  fluid  dynamics)  when  only  moderate 
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accuracy  was  required  [1].  More  recently,  an  integrator  based  heavily  on  CHEMEQ  outperformed  a  first-order 
quasi-steady-state  method  and  the  implicit  preconditioning  method  CHEMSODE  [30]  on  a  photochemical 
smog  problem  [29].  Despite  its  strengths,  CHEMEQ  exhibits  instability  under  some  situations  and  is  limited 
in  the  accuracy  it  can  achieve  [19]. 

This  report  describes  a  quasi-steady-state  method  which  we  call  a-QSS,  and  its  implementation  in  the 
subroutine  CHEMEQ2.  The  a-QSS  method  is  A-stable  for  linear  problems  and  second-order  accurate.  It 
is  more  stable,  more  accurate,  and  costs  less  than  CHEMEQ,  and  it  successfully  integrates  some  systems 
for  which  CHEMEQ  fails  [19],  CHEMEQ2,  has  been  used  successfully  in  hydrogen- air  flame  studies  in  mi¬ 
crogravity  [16],  on  pulse-detonation  engine  studies  [17, 18],  on  thermonuclear  mechanisms  used  in  supernova 
simulations,  and  on  the  test  cases  used  to  validate  CHEMEQ  [19],  In  addition  to  describing  the  new  algo¬ 
rithm,  we  present  error  and  linear  stability  analyses.  We  also  describe  how  to  use  the  subroutine  CHEMEQ2 
as  a  stand-alone  integrator  and  in  conjunction  with  a  reacting-flow  code.  Variables  used  in  CHEMEQ2 
are  listed  and  documented  in  Appendix  A,  and  results  obtained  using  CHEMEQ2  are  compared  to  those 
obtained  using  CHEMEQ  on  two  test  problems.  Finally,  code  listings  for  CHEMEQ2  and  its  supporting 
subroutines  are  also  included, 

2  Introduction  to  QSS  Methods 

Consider  a  simplified  form  of  Eq.  (1),  in  which  the  subscript  i  is  dropped  for  convenience,  =  0,  and 

y(<”)  =  y\ 

^  =  q-py  y(0)  =  y‘’.  (3) 

If  p  and  q  are  constant,  then  Eq.  (3)  has  an  exact  solution  given  by 

y(<)  =  ^  (l  -  e"'’*) .  (4) 

Quasi- steady- state  (QSS)  methods  are  based  on  the  solution  given  in  Eq.  (4)  [21-24].  If  q  and  p  are  slowly 
varying,  evaluating  Eq.  (4)  a,t  t  =  At  using  q{i^)  and  p{t^)  provides  a  good  approximation  for  y{At).  This 
approach  gives  a  first-order  method  which  is  the  simplest  QSS  algorithm.  More  sophisticated  QSS  algorithms 
incorporate  the  time  dependence  of  p  and  q  and  may  place  Eq.  (4)  into  an  alternate  algebraic  form.  The 
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common  thread  between  QSS  methods  is  their  basis  on  Eq.  (4),  which  requires  the  methods  to  return  the 
exact  solution  if  9  and  p  are  constant.  There  are  many  QSS  methods  documented  in  the  literature,  and  the 
or-QSS  method  is  compared  to  several  of  them  in  Section  3.2. 

Note  that  if  pAt  00  (i.e.,  the  ODE  timescale  is  very  small  compared  to  the  numerical  timestep).  the 
exponential  terms  in  Eq.  (4)  may  be  neglected  completely,  leading  to  the  steady-siate  approximation  [31] 


Since  g  and  p  are  functions  oft,  a  steady-state  approximation  for  species  f  does  not  imply  that  the  value 
of  y.  remains  constant.  Equation  (5)  assumes  that  the  adjustment  toward  a  “local”  equilibrium  for  species 
t  based  on  the  current  values  of  g,-  and  p,-  is  instantaneous.  In  contrast,  quasi-steady-state  methods  retain 
information  about  the  transition  to  equilibrium  and  are  therefore  more  accurate  for  larger  timescales.  Note 
that  some  authors  call  Eq.  (5)  a  quasi-steady-state  approximation  [32],  emphasizing  the  continued  evolution 

of  y.-  as  g,  and  p,  change.  In  this  paper,  the  label  “quasi-steady  state”  is  reserved  for  methods  that  reproduce 
the  solution  in  Eq.  (4)  for  constant  g  and  p  regardless  of  the  timescale. 

QSS  methods  are  often  compared  with  standard  stiff  solvers  such  as  LSODE  [25,26],  which  is  a  variable- 
order  method  based  on  Gear’s  backward  differentiation  formulae  (BDF)  [27].  However,  such  comparisons 
have  been  largely  limited  to  the  integration  of  a  single  problem  from  one  set  of  initial  conditions,  not 
reacting-flow  simulations  in  which  start-up  overhead  and  storage  requirements  play  key  roles  in  the  overall 
efficiency  of  the  integrator.  Verwer  and  Simpson  describe  one  such  test  from  atmospheric  chemistry,  in  which 
a  simple  twcvstep  BDF  method  outperforms  a  first-order  implicit  QSS  method  and  a  two-stage  explicit  QSS 
method.  The  test  involved  the  calculation  of  emissions  and  was  not  coupled  to  fluid  dynamics  [22].  Jay  et  al. 
introduce  two  QSS  methods  and  examine  their  performance  on  a  set  of  atmospheric  tests  involving  32  species 
[21].  These  two  QSS  methods  outperformed  both  a  standard,  first-order  QSS  method  and  CHEMEQ,  but 
the  methods  were  slower  than  multi-point  BDF  methods.  Variable-order,  multi-point  BDF  methods  often 
outperform  QSS  methods  when  the  chemistry  integration  stands  alone.  However,  the  demands  of  a  reacting- 
flow  application  are  very  different  than  tho.se  of  a  stand-alone  integration,  and  the  conclusions  of  these  studies 
cannot  be  applied  to  reacting-flow  problems.  The  o-QSS  algorithm  was  developed  specifically  to  meet  the 
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demands  of  a  process-split,  reacting-flow  simulation. 


3  The  a-QSS  Algorithm 

3.1  Algorithm  Development 


Given  the  demands  of  a  reacting-flow  application,  we  chose  a  predictor-corrector  implementation  for  the 
integrator.  Evaluating  Eq.  (4)  at  At  using  initial  values  serves  as  the  predictor  step,  and  a  correction  based 
on  the  initial  and  the  predicted  values  then  follows.  The  corrector  step  can  be  repeated  using  the  previous 
corrector  result  as  the  new  predicted  value.  Predictor-corrector  methods  of  this  type  are  single-point  methods 
because  information  from  only  a  single  time  level  is  needed  to  initiate  calculation  of  the  solution  at  the  next 
time  level. 

First,  a  convenient  algebraic  form  for  Eq.  (4)  was  chosen.  Equation  (4)  can  be  evaluated  a,t  t  —  At, 
yielding 


y{At)  =  y°  + 


At{q-py°) 

1  -|-  ocpAt  ’ 


(6) 


for  a  defined  by 


a(p  At)  = 


l_(l_e-pA>)/(pAt) 
1  -  e-f  A< 


(7) 


The  parameter  a  is  a  function  of  p  At,  as  shown  in  Fig.  1.  Note  that  a  0  as  p  At  — -oo,  a  — ►  1  as 
p  At  oo,  and  a  =  1/2  for  p  At  =  0.  The  meanings  of  these  limits  are  clarified  by  recalling  that  p  At  —  At /r. 
The  a  ^  1  limit  corresponds  to  an  infinitely  fast  ODE  relative  to  At,  and  a  =  1/2  corresponds  to  an  infinitely 
slow  ODE.  Equation  (6)  is  exact  for  any  value  of  p  (provided  q  and  p  are  constant).  However,  we  split  g  such 
that  pp  is  a  non-negative  loss  rate,  so  only  values  of  p  At  >  0  need  be  considered.  Approximations  used  to 
calculate  a(pAt)  without  the  costly  exponential  evaluation  are  described  in  Section  4.1, 

A  predictor-corrector  method  based  on  the  solution  in  Eq.  (6)  takes  the  form 


if  ~y^  ^ 


At  (g^-pV) 
1-f  oO  At  po 


Predictor, 


(8) 


At  {q*-p*y*) 
1  a  *  At  p* 


Corrector. 


(9) 


5 


-50 


0 

p  At 


50 


100 


Figure  1;  The  parameter  a  as  a  function  of  pAt. 
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Superscript  0  indicates  initial  values,  and  superscripts  p  and  c  indicate  predicted  and  corrected  values, 
respectively.  The  predictor  uses  the  initial  values  of  q,  p,  and  y,  but  the  ''starred”  variables  (g*,  p*,  y*,  and 
a*)  can  be  based  on  both  the  initial  values  and  the  predicted  values. 

If  we  assume  linear  profiles  in  time  for  q  and  p  between  the  initial  and  predicted  values,  we  can  find 
an  exact  series  solution  for  Eq.  (1).  (This  solution  is  illustrated  in  conjunction  with  the  error  analysis  in 
Section  3.3.)  Unfortunately,  the  series  solution  does  not  readily  provide  an  efficient  integration  technique,  nor 
does  it  indicate  appropriate  averages  for  the  starred  variables  in  the  corrector.  However,  solutions  do  exist 
under  slightly  simpler  conditions  that  can  be  reproduced  with  appropriate  choices  of  the  starred  variables. 

For  instance,  if  p  is  constant  and  q  is  linear  in  time,  the  exact  solution  to  Eq.  (1)  can  be  written  as 


y(A<)  =  y®  + 


At  {q-  py^) 

1  +  a  At  p  ’ 


(10) 


for  a  =  a{pAt)  from  Eq.  (7)  and 


q  =  aq{At)  +  {1  -  a)q^ 


Alternatively,  if  5  =  0  and  p  is  linear  in  time,  the  exact  solution  of  Eq.  (1)  is 


in  which 


(11) 


(12) 


P  =  ^  (p(A<)+p")  , 


(1.3) 


and  a  =  a(pAt)  from  Eq.  (7). 

These  results  suggest  a  corrector  of  the  form 


y"  + 


{q-py°) 

l-f-oAtp 


(14) 


To  calculate  q  and  p  from  Eqs.  (11)  and  (13),  we  replace  ^(At)  and  p(At)  with  the  predicted  values  q^  and 
pP.  When  q  and  p  are  known  functions  of  the  predicted  values  are  replaced  with  the  exact  values  at  At. 
Using  {^,p,y^,a}  for  y*, a"')  in  Eq.  (9)  gives  a  method  which  is  A-stable  for  linear  problems  and 

second-order  accurate,  as  is  shown  in  Sections  3.3  and  3.4.  We  refer  to  the  new  method  as  a-QSS,  which 
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emphasizes  the  dual  role  that  a  plays  in  returning  the  exact  solution  for  constant  q  and  p  and  in  providing 
a  weighted  average  of  q  when  q  is  not  constant. 

3.2  Comparison  to  Previous  Methods 

In  addition  to  the  algebraic  form  chosen  for  Eqs.  (8)  and  (14),  a-QSS  differs  from  previous  QSS  methods 
in  its  choice  of  averaging  and  its  implementation  as  a  predictor-corrector  method.  Previous  methods  that 
calculate  average  values  for  p  and  q  use  the  same  averaging  method  for  both  terms.  For  example,  the  two- 
stage  explicit  method  introduced  by  Verwer  and  Van  Loon  [23]  and  tested  by  Verwer  and  Simpson  [22]  uses 
a  simple  algebraic  average  for  both  q  and  p  calculated  from  initial  and  predicted  values.  CREKID  [24]  uses 

an  implicit  exponential  Euler  formulation  in  which  a(pAt)  gives  a  weighted  average  of  the  composite  source 
terms; 


y(At)  —  t/*  +  At  (ag(At)  +  (1  —  ar)<jr*’)  .  (15^ 

In  contrast,  the  a-QSS  algorithm  uses  a  simple  algebraic  average  for  p  and  an  o-weighted  average  for  q  in 
order  to  match  the  exact  solutions  described  in  Eqs.  (10)  and  (12). 

Other  QSS  methods  combine  the  results  of  first-order  calculations  in  a  way  that  improves  accuracy.  Jay 
et.  al.  [21]  describes  two  such  methods.  Their  “extrapolated  QSS”  method  finds  the  solution  at  t°+At,  first 
with  a  single  step  and  then  with  two  steps  of  At/2  each.  A  simple  extrapolation  then  estimates  the  solution 
that  would  result  if  an  infinitely  small  timestep  were  used.  Their  second  method,  “symmetric  QSS,”  is  a 
two-step  method  requiring  three  evaluations  of  the  source  terms.  Each  of  these  steps  acts  as  if  9  and  p  were 
constant,  and  the  values  for  q  and  p  are  taken  at  the  same  time  level  based  on  the  previous  calculation.  No 
averaging  of  ^  or  of  p  occurs  between  time  levels  in  these  methods. 

The  algebraic  form  of  Eqs.  (8)  and  (14),  which  was  introduced  in  Eq.  (6),  is  based  on  the  asymptotic 
update  employed  by  CHEMEQ  when  the  timescale  for  an  equation  is  smaller  than  some  user-specified  value 

[1,2].  However,  CHEMEQ  effectively  replaces  o(pA<)  with  the  constant  1/2,  which  is  equivalent  to  choosing 
the  Fade  approximation 


exp(a;) 


2  +  X 
2- X 


(16) 
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in  either  the  definition  of  a  or  in  Eq.  (4).  When  the  timescale  for  an  equation  is  larger  than  some  user-specified 
value,  that  equation  is  integrated  using  the  modified  Euler  method.  The  hybrid  method  studied  by  Lorenzini 
and  Passoni  [29]  uses  CHEMEQ’s  update  equations  but  different  criteria  for  determining  the  timestep  and 
for  choosing  between  the  asymptotic  update  and  the  modified  Euler  update.  CHEMEQ’s  asymptotic  update 
also  uses  different  averages  in  the  corrector  for  p  and  q  than  those  used  in  Eq.  (14).  These  differences  lead 
to  instability  in  CHEMEQ  that  is  discussed  more  thoroughly  in  Section  6.1,  The  averages  chosen  by  or-QSS 
eliminate  this  instability,  and  a-QSS  automatically  approaches  the  modified  Euler  method  as  pAt  0. 

3.3  Error  Analysis 


The  method  has  a  third-order  error  term  for  a  single  step,  which  makes  it  second-order  over  the  course  of 
an  integration.  This  can  be  shown  by  examining  the  exact  series  solution  of  Eq.  (1),  Writing  the  series  for 
y{t)  about  y{t^  =  0)  =  t/Oi 


the  derivative  is  given  by 


y(t)  =  yo  +  yit  + +  •  •  •  = 

3=0 


^  =  yi  +  2y2t +  +Zy3t^  +  ■■■  =  ^jyji^  ^ 

j  =  l 


(17) 


(18) 


This  development  deals  with  a  single  species,  j/,  so  yj  is  the  coefficient  of  the  P  term  in  the  expansion  in 
Eq.  (17)  and  not  the  concentration  of  the  species  in  a  multi-species  system.  Similarly,  series  expansions 
for  q{t)  and  p(<)  are  given  by 


9(0  = 


3=0 


p(0  =  X!pjO. 

3=0 


(19) 


(20) 


Substitution  into  Eq.  (1)  using  Eqs.  (17)-(20)  gives 


yi  =  qo  -  poyo, 


1/2  =  2  (91  -  (PiVo  +  PoVi))  , 

ya  =  ^  (92  -  (i>2yo  +  piyi  +  poy^)) , 


(21) 

(22) 

(2.3) 
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and  leads  to  the  general  expression 


for  j  >  0.  (The  exact  solution  for  q  and  p  both  linear  in  time  that  was  mentioned  in  Section  3.1  is  obtained 
by  assuming  pj  -  qj  -  0  for  j  >  2.) 

In  general,  q  and  p  are  given  as  functions  of  y,  not  as  functions  of  t.  Therefore,  the  coefficients  in 
Eqs.  (19)  and  (20)  are  not  known,  and  Eq.  (3)  is  a  nonlinear  differential  equation.  We  will  first  perform  an 
error  analysis  for  the  linear  version  of  Eq.  (3),  in  which  q  and  p  are  known  functions  of  t,  and  then  extend 
this  analysis  for  the  nonlinear  case.  For  the  linear  case,  the  predicted  values  are  simply  qP  =  q{At)  and 
pP  =  p{At).  Subtracting  the  series  expansion  for  Eq.  (14)  from  the  exact  solution  evaluated  at  t  =  At  yields 

y{At)  -  y'  =  _  52  +P2!/o^ 

•f  0(A^^)  [linear  case]  ^25) 

The  leading  error  term  is  0(At3)  per  timestep.  Since  the  number  of  timesteps  required  to  reach  a  given  time 
is  proportional  to  1/At,  the  error  for  the  method  is  second-order  when  these  errors  all  reinforce  [3].  Note 
that  this  error  term  disappears  for  constant  p  and  linear  5  (which  gives  pi  =  52  =  p2  =  0)  and  linear  p  with 

q  =  0  (which  gives  50  =  92  =  P2  =  0).  These  two  cases  were  used  to  choose  how  to  calculate  the  starred 
variables  in  Eq.  (9),  and  the  method  is  exact  for  either  case. 

The  method  is  second-order  for  nonlinear  problems  as  well.  To  illustrate  this,  first  note  that  the  leading 
error  term  for  the  predicted  values  is  second-order: 


y{At)  -  yP  =  -^(91  -  piyo)  +  O(At^).  (26) 

Since  9  and  p  are  polynomials  in  the  species  concentrations  for  the  nonlinear  systems  representing  reaction 

kinetics,  the  leading  error  terms  for  the  predicted  values  qP  and  pP  are  also  second-order.  This  error  can  be 
represented  as 


q{At)  -  qP  =  (gAt-  +  0{At^), 


(27) 


p(A<)  -  ;/  =  epAt-  +  0(At3), 


(28) 
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for  some  unknown  coefficients  €q  and  Cp.  Using  these  predicted  values  in  Eq.  (14)  gives  an  error  term  of  the 
form 


yiAt)-y^  =  At^  (-1 

-h  0{At'^)  [nonlinear  case]. 


11  11 
PiQo  g^2  4-  gP2yo  4-  —  ^^pVO 


(29) 


As  with  the  linear  problem,  the  leading-order  error  term  for  the  nonlinear  problem  is  0(A/^)  per  timestep, 
so  the  method  is  still  at  least  second-order  over  the  course  of  an  integration. 

3.4  Linear  Stability  Analysis 


For  the  single  linear  equation 


dt 


the  coefficient  A  can  be  a  function  of  t  but  not  a  function  of  y.  Using  the  average  value  A  given  by 

A  =  i(A(t  =  0)  +  A(t  =  A<)), 


a-QSS  has  amplification  factor  G  given  by 


G=  1  + 


XAt 


(30) 


(31) 


(32) 


1  —  a  A  A/ 

The  signs  in  Eq.  (32)  reflect  the  fact  that  A  =  — p,  and  note  that  a  =  a(—XAt).  Using  Eq.  (7),  the  expression 
for  G  simplifies  to 


G  =  exp(A  At). 


For  A  =  a  -f  6\/— 1  with  a,  6  both  real,  the  magnitude  of  G  is  simply 


(33) 


||G||  =  exp  {a  At) . 


(34) 


Since  |1G||  <  1  for  a  <  0  for  any  value  of  6,  the  method  is  A-stable.  This  does  not  prove  that  q-QSS  is 
A-stable  when  applied  to  nonlinear  systems  of  ODEs  for  which  {pi}  and  {qi]  depend  on  {pt}*  However, 
in  testing  to  date,  the  accuracy-based  timestep  criterion  used  originally  in  CHEMEQ  has  worked  well  for 
o-QSS.  To  ensure  that  stability  is  maintained  when  the  corrector  is  iterated,  a  new  stability  check  was 
introduced  [20]  for  A^.  These  accuracy  and  stability  criteria  are  discussed  in  Section  4.2. 
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4  CHEMEQ2  Implementation 

4.1  Update  Equations 

CHEMEQ2  uses  the  a-QSS  update  on  all  equations  in  the  system  regardless  of  the  timescale  of  the  ODE. 

This  makes  the  timestep  less  prone  to  oscillations  than  for  hybrid  methods  (like  CHEMEQ)  that  switch 

between  update  methods.  In  addition,  iterations  may  be  done  on  the  corrector  that  improve  the  accuracy 
of  the  result. 

Again  using  a  superscript  0  to  indicate  values  at  the  begining  of  the  chemical  timestep  and  a  subscript  i 
to  specify  species  t,  the  a-QSS  update  is  given  by 


/  =  + 


Predictor, 


(35) 


Vi 


0  , 

l+^AtpI 


Corrector. 


(36) 


The  predictor  uses  all  initial  values,  and  a?  =  a(p?A0.  After  calculating  the  predicted  concentrations  {yf } 
for  all  of  the  species  in  the  system,  next  obtain  {gf }  and  {pf }  from  {y^}.  Then  calculate 


P*  —  2  > 

evaluate  57  =  a(^At),  and  finally 


(37) 


ft  -  0!iq^  +  (1  -  a,)g9 . 

These  averages  are  then  used  to  calculate  yf,  and  to  iterate  on  the  corrector,  use  the  value  for  yf  from  one 
step  as  for  the  next  step. 

Having  an  accurate  approximation  for  a(pAt)  that  does  not  require  an  evaluation  of  the  exponential 
function  makes  the  method  given  by  Eqs.  (35)  and  (36)  more  attractive.  Recall  that  p  is  strictly  non-negative 
based  on  the  way  the  chemical  source  term  is  split,  so  this  approximation  need  only  hold  for  positive  values 
of  p  At.  Equation  (7)  indicates  that  as  pAt  — ^  oo,  a  reasonable  approximation  for  a(pAt)  is 


Qf(p  A/.)  I - 

pAt 


(39) 
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Using  this  approximation  for  a  eliminates  the  need  to  find  an  accurate  approximation  for  as  x  — j*  oo,  as 

would  be  required  if  the  solver  were  based  on  Eq.  (4)  rather  than  Eq.  (6),  Using  the  Fade  approximation  [4] 

^  ^  1680  +  840x  +  180x2  +  20x3  +  x^ 

^  ~  1680  -  840a;  +  180a;2  -  20a;3  + 


in  the  definition  of  a{pAt)  gives 


a(Af/r) 


840r3  +  140r2  +  20r  +  1 
1680r3  +  40r 


(41) 


for  r  =  l/(pAf).  These  two  approximations  are  shown  with  the  exact  curve  for  a  in  Fig.  2.  The  approxi¬ 
mation  given  by  Eq.  (41)  is  labeled  Fade  (a).  Note  that  unlike  Fig.  1,  tbe  x-coordinate  in  Fig.  2  is  r/Ai. 
The  linear  approximation  in  Eq.  (39)  is  closer  to  the  exact  value  of  a  than  the  approximation  in  Eq.  (41)  for 
r  <  0.16762;  for  r  >  0.16762,  Eq.  (41)  is  more  accurate.  Therefore,  the  better  approximation  can  be  chosen 
based  upon  the  value  of  r.  This  combined  approximation  differs  from  the  exact  value  of  a  by  at  most  0.3%. 
This  error  occurs  in  a  narrow  region  around  the  transition  from  the  linear  to  the  rational  approximation. 
As  r  — oo,  Eq.  (41)  is  actually  better  behaved  numerically  than  Eq.  (7),  which  requires  the  exponential 
function  evaluation. 
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A  second  approximation  for  a  that  can  be  used  alone  (i.e.,  does  not  require  the  linear  approximation 
Eq.  (39),  to  take  over  as  r  ^  0)  results  from  the  Fade  approximation 


^ _ ^0  +  120a;  +  12x2 

~  360  -  240x  +  72x2  -  12x3  +  ' 


(42) 


This  approximation  gives 


a(A</r) 


180r3  +  60r2+llr+l 


360r3  +  GOr^  +  12r  +  1 '  (^3) 

This  second  approximation  for  a  is  designated  Fade  (b)  and  is  also  shown  in  Fig.  2.  Since  this  approximation 
recovers  a  =  1  for  r/A<  =  0,  Eq.  (43)  can  be  used  alone  with  only  a  slight  accuracy  penalty  compared  to 
the  combined  approximation  of  Eqs.  (39)  and  (41).  In  testing  to  date,  this  accuracy  penalty  in  a  has  not 
caused  accuracy  problems  in  the  species  solutions,  and  using  this  single  equation  eliminates  the  logic  check 
required  to  determine  which  of  Eqs.  (39)  or  (41)  to  use  for  the  combined  method.  The  success  of  CHEMEQ, 
which  effectively  uses  a  =  1/2,  suggests  that  an  even  simpler  approximation  for  a{pAt)  than  those  shown 
here  may  provide  sufficient  accuracy  with  lower  computational  cost. 


4.2  Timestep  Selection 

Accuracy  is  controlled  by  choosing  At  and  the  number  of  corrector  iterations,  N,.  A  single  corrector 
calculation  is  performed  if  AT,  =  1,  and  as  Appendix  B  illustrates,  increasing  N,  improves  accuracy.  Timestep 

selection  is  identical  to  that  used  by  CHEMEQ  [2].  The  initial  predicted  values  and  final  corrected  values 
are  tested  to  see  if  they  satisfy 


!l2/i  (44) 

for  some  specified  constant  c,  typically  ~  10-2.  If  Eq.  (44)  is  not  satisfied,  the  step  is  repeated  with  a  smaller 
timestep.  As  Young  notes  [2],  it  is  best  to  reduce  the  timestep  sharply  (a  factor  of  2  or  3)  rather  than  slowly 
as  less  time  is  wasted  finding  a  step  size  sufficiently  small  for  convergence.  However,  if  the  convergence 
criterion  is  met  easily  during  the  iteration,  it  is  best  to  only  increase  the  step  size  by  5-10%  each  step  [2]. 
The  timestep  modification  is  performed  by  modeling  the  difference  between  predictor  and  corrector  as  a 
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single  second-order  term.  Choose  a2  such  that 


Vi  -  Vi  =  a2(A<)L.  (45) 

where  (A^)om  is  the  timestep  used  to  calculate  and  yj  from  the  initial  conditions.  The  user  specifies  a 
target  value  for  the  relative  magnitude  of  this  correction  term,  given  by 

||a2A<^||tarje<  =  £!/,••  (46) 

The  values  of  e  in  Eq.  (44)  and  e  in  Eq.  (46)  are  related  by 

€  =  C£  (47) 


for  some  user-specified  c  >  1,  so  the  error  criterion  in  Eq.  (44)  is  rarely  violated  for  the  chosen  timestep. 
Defining  the  parameter  cr  as 


cr  =  max 

a 


( 


eyf  /  ’ 


(48) 


the  value  of  At  that  limits  the  largest  relative  concentration  change  to  the  target  magnitude  is 

(AOtar^et  =  (49) 

y/o- 

Again,  the  difference  \\yf  —  yf  ||  is  calculated  using  the  initial  prediction  for  y^  and  the  final  corrected  value 
oft/f. 


In  CHEMEQ2,  the  new  timestep  is  calculated  using 


(A<)ne«.  -  {At)old  +  .005 j  .  (50) 

In  Eq.  (50),  is  an  estimate  by  three  Newton  iterations  for  y/a  using  a  as  the  starting  value.  Equation  (50) 
gives  the  desired  asymmetrical  property  in  that  At  decreases  faster  than  At  would  increase  for  the  inverse 
value  of  O'.  In  addition,  At  is  modified  very  little  when  a  is  near  1, 

If  Eq.  (44)  is  satisfied,  the  new  concentration  values  at  i  -f  (A^)o/d  are  set  equal  to  the  values  of 
y^  and  the  integration  continues  with  timestep  (A/)n€ty*  A  successful  integration  step  requires  only  two 
derivative  function  evaluations  when  a  single  corrector  step  is  used,  and  the  timestep  selection  algorithm 
minimizes  the  number  of  steps  that  must  be  repeated. 
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When  yi  is  decaying  toward  zero,  it  is  constrained  by  a  minimum  value.  When  this  lower  bound  is 
reached,  the  species  is  no  longer  considered  in  the  calculation  of  a  and  therefore  does  not  affect  convergence. 
The  lower  bound  is  chosen  by  the  user  and  must  be  selected  to  insure  that  it  does  not  adversely  affect  the 
accuracy  of  the  solution.  If  the  minimum  value  is  too  high,  a  species  value  could  freeze  prematurely.  This 
could  corrupt  the  solution  for  other  species  whose  values  are  sensitive  to  the  frozen  species.  If  the  minimum 
is  too  low,  it  can  affect  efficiency  if  the  decaying  mode  controls  the  timestep.  Thus  it  is  better  to  choose 
minimum  values  that  are  too  small  rather  than  too  big.  This  may  slightly  reduce  computational  efficiency, 
but  It  will  also  reduce  the  possibility  of  calculating  an  incorrect  solution. 

As  mentioned  in  Section  3.4,  a-QSS  is  A-stable  for  linear  problems,  but  this  result  holds  no  guarantees 
for  the  nonlinear  systems  of  chemical  kinetics.  To  ensure  that  the  calculation  remains  stable,  the  integrator 
can  monitor  the  convergence  of  the  corrector  iteration  and  adjust  the  timestep  if  necessary.  Let  yf^  denote 
the  corrected  value  of  y.(AO  after  /  iterations.  The  change  from  one  iteration  to  the  next, 

CO-1) 

should  decrease  in  size  as  /  grows  if  the  iteration  is  stable.  Therefore,  requiring  that 


||Aj,r‘1|<||Ay‘ 


(52) 


where  N,  is  the  specified  maximum  number  of  iterations,  ensures  that  the  integration  remains  stable.  A  step 

using  At  that  satisfies  the  accuracy  constraint  but  fails  to  satisfy  Eq.  (52)  for  any  i  is  repeated  using  a  new 
timestep,  given  by 

(53) 


(At)„e«,  =  At  max  ( — 

‘  \l|Aj/f^'^||  + 0.001 


This  instability  is  generally  not  seen  in  reacting-flow  applications,  so  a  more  sophisticated  criterion  and 

timestep  update  have  not  been  pursued.  The  constraint  is  available  in  CHEMEQ2,  however,  should  it  be 
needed- 

At  start-up,  an  initial  timestep  is  chosen  which  approximates  the  timestep  that  will  be  determined  by 
the  predictor-corrector  scheme.  This  initial  trial  timestep  is  determined  sucli  that  none  of  the  variables  will 
change  by  more  than  a  prescribed  amount.  If  the  formation  rate  ,,  is  much  larger  than  the  loss  rate  p.  y,. 


16 


it  is  reasonable  to  assume  that  qi  and  p,  will  remain  relatively  constant  for  large  changes  in  y,- .  The  initial 
change  in  y,  may  be  large  enough  to  equilibrate  the  formation  and  loss  rates  for  yi ,  Thus  the  initial  trial 
timestep  is  chosen  as 

At  =  e  min  or  (if  qi  >  pjt/i)  1/pi^  ,  (54) 

where  e  is  a  scale  factor  that  need  not  be  identical  to  the  constant  used  in  Eq.  (44).  Equation  (54)  is  used  only 
once  per  global  timestep,  as  subsequent  timesteps  taken  until  the  end  of  the  global  timestep  are  determined 
using  Eq.  (50). 

The  effect  of  the  thermodynamic  state  on  the  the  reaction  rate  constants  has  been  ignored  in  the  previous 
developments.  The  rate  constants  are  often  calculated  once  before  the  predictor  step  using  the  initial  values 
and  held  constant  during  the  corrector  step.  A  new  thermodynamic  state  is  then  found  at  the  new  time 
level  and  used  for  the  following  predictor  and  corrector.  If  the  integration  is  particularly  sensitive  to  the 
thermodynamic  state,  this  state  can  be  recalculated  for  the  corrector  using  the  predicted  solution.  If  the 
system  requires  the  integration  of  a  thermodynamic  variable  (such  as  temperature)  along  with  the  species 
concentrations,  then  the  source  term  for  this  extra  variable  is  split  just  as  with  the  concentrations.  If  there  is 
no  ‘^loss’’  term  for  that  variable  that  can  be  assumed  proportional  to  the  variable,  then  the  entire  source  term 
is  assigned  to  y,  and  the  method  reduces  to  the  modified  Euler  method  for  that  equation  since  o;(0)  =  1/2. 

5  How  to  Use  CHEMEQ2 

The  following  sections  describe  what  a  user  must  know  about  CHEMEQ2  in  order  to  use  the  subroutine 
effectively.  A  description  of  the  four  entrance  points  for  the  subroutine  and  the  argument  variables  used 
for  each  is  included,  as  are  diagrams  indicating  the  calling  sequence  used  in  a  stand-alone  integration  and 
in  a  reacting-flow  application.  Appendix  A  describes  all  internal  variables,  and  Appendix  B  provides  code 
listings, 

CHEMEQ2  is  designed  as  a  replacement  for  CHEMEQ,  so  most  of  the  original  CHEMEQ  code  is  retained 
in  CHEMEQ2.  The  overall  logical  structure  from  CHEMEQ  is  retained,  with  the  hybrid  method  replaced 
with  a-QSS-  Minor  changes  in  input/output  have  also  been  implemented,  and  these  are  discussed  in  the 
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following  section. 


5.1  Entry  Points 

1.  CHEMEq2(<ltg.gsub.n.y)  is  the  main  entrance  point  and  is  used  to  advance  the  chemical  variables  in 
time. 


•  dtg:  time  for  which  the  integration  is  performed;  Atg  of  Section  1. 

•  gsub:  name  of  the  derivative  function  evaluator  subroutine  that  provides  the  source  term  as  q 
and  py.  The  form  of  gsub  and  its  arguments  are  given  in  Section  5.2. 

•  n:  then  number  of  equations  to  be  integrated 

•  y:  array  which  holds  the  initial  values  at  the  start  of  the  integration  and  returns  the  final  values 
at  the  end  of  the  integration 

2.  CHEMSPfepsmn.  epsMx.  dtnui.  tuot.  iterrux.  us.  yam.  prt)  provides  a  means  to  set  the  solution 
parameters  used  the  next  time  CHEMEQ  is  called.  If  the  passed  variable  has  value  <  0,  then  the 
default  value  built  into  the  subroutine  is  used.  If  the  passed  value  is  >  0,  then  the  corresponding 
parameter  is  set  to  the  passed  value.  For  a  typical  calculation  in  which  the  same  solution  parameters 
are  used  throughout  the  domain,  this  routine  may  be  called  only  once  to  initialize  these  parameters.  If 
the  simulation  involves  multiple  regions  that  make  different  demands  on  the  speed  or  accuracy  of  the 
integration,  then  this  routine  may  be  called  so  that  appropriate  parameters  are  used  in  each  region. 
The  parameters  set  by  each  variable  and  the  default  value  for  these  parameters  are  listed  below.  The 
distinction  is  made  between  the  arguments  of  CHEMSP  and  the  internal  variables. 

.  epsmn:  sets  epsmin,  the  accuracy  parameter  s  in  Eq.  (46)  for  determining  the  timestep.  Default 
value  of  epsmin:  10“^. 

.  epsmx:  sets  epsmax,  the  parameter  c  in  Eq.  (47)  for  specifying  when  a  step  must  be  repeated 
using  a  smaller  timestep.  Default  value  of  epsmax:  10. 

•  dtmn:  sets  dtmin,  the  minimum  allowed  timestep.  Default  value  of  dtmin:  IQ-*®. 
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•  tnot:  sets  t start,  the  value  of  the  independent  variable  at  the  start  of  the  global  timestep. 

Default  value  of  t start:  0. 

•  itermx:  sets  itermaoc,  the  number  of  corrector  iterations  performed.  Default  value  of  itermajc: 

1  (a  single  corrector  step). 

•  ns:  integer  that  indicates  the  number  of  entries  in  ymin  to  initialize  with  yum 

•  ynm:  sets  ymin,  the  array  which  holds  the  minimum  allowed  values  of  the  dependent  variables  (for 

integration  control).  Default  value  of  ymin(i):  10“^®  for  all  values  of  i. 

•  prt:  if  prt  =  0,  the  list  of  parameters  set  by  the  current  call  to  CHEMSP  is  printed. 

3.  CHEMCT(tink)  provides  diagnostics  by  printing  the  number  of  derivative  function  calls  and  the  number 
of  times  an  integration  step  was  redone  due  to  a  violation  of  the  accuracy  criterion  or  the  stability 
criterion, 

•  tmk:  REAL  number  used  to  identify  the  call;  the  value  of  the  independent  variable  is  often  used. 

4.  CHEMPR(y  ,n)  is  called  for  diagnostic  purposes  and  prints  partial  lists  of  internal  variables.  The  variable 
definitions  are  the  same  as  those  in  the  CHEMEQ2  call. 

5.2  Supporting  Subroutines 

1.  CHEMER  is  a  diagnostic  routine  which  warns  the  user  that  the  minimum  timestep  threshhold  has 
been  violated  and  that  the  integration  has  been  stopped.  No  arguments  are  required.  The  user  may 
supply  additional  error-checking  capabilities  or  diagnostic  output. 

2.  gsub(y,  q,  d,  t)  is  the  derivative  function  evaluator.  The  actual  name  of  this  subroutine  is  passed 
to  CHEMEQ2  as  an  argument.  This  is  so  the  routine  can  be  changed  in  different  regions  or  regimes. 

•  y:  array  holding  the  dependent  variables 

•  q:  production  terms;  entry  q(i)=  in  Eq.  (1). 

•  d:  loss  terms;  entry  d{±)=  piyi  in  Eq.  (1). 

•  t:  current  value  of  the  independent  variable 
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Figure  3;  Basic  calling  sequence  to  obtain  final  values 


5.3  Calling  Sequence 

If  the  user  simply  needs  the  species  concentrations  at  for  a  single  problem,  the  basic  calling  sequence 
illustrated  in  Fig.  3  should  be  used.  A  single  call  to  CHEMSP  is  required  to  initialize  the  integration 
parameters  unless  all  default  values  are  used,  and  dtg  =  Using  the  version  of  CHEMEQ2  included  in 

this  report,  no  intermediate  values  are  provided  between  t  =  0  and  The  driver  program  included  in 

Appendix  C  has  this  structure  imbedded  in  a  loop  that  scrolls  through  various  solution  parameters  to  give 
a  set  of  final  values.  If  intermediate  values  are  desired,  a  write  statement  can  be  added  within  CHEMEQ2 
in  order  to  print  the  results  of  a  successful  step  before  the  next  step  is  taken.  To  obtain  intermediate  values 
without  altering  CHEMEQ2,  the  integration  time  tjinai  may  be  broken-up  as  illustrated  in  Fig.  4.  After  each 
At,  step,  control  returns  to  the  driver  program  and  values  may  be  printed.  The  next  step  in  the  integration 
is  taken,  and  the  result  of  the  previous  step  is  used  as  the  initial  condition.  The  optional  CHEMSP  call  is 
needed  if  the  source  terms  depend  on  i  (such  as  in  atmospheric  chemistry)  or  if  the  solution  parameters  are 
to  be  changed  (i.e.,  perhaps  the  initial  time  interval  is  very  sensitive  and  must  be  run  more  accurately,  but 
later  times  allow  this  constraint  to  be  loosened).  Variables  changed  via  CHEMSP  and  passed  as  arguments 
to  CHEMEQ2  may  be  moved  if  necessary  for  efficiency.  Please  see  the  Section  5.4  for  a  discussion  of  some 
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Calculate  initial  conditions, 
control  parameters 

1 

Call  CHEMSP 

1 

setAt^ 

; - -1 

optional  CALL  CHEMSP 


Call  CHEMEQ2 


Print  diagnostics. 


Figure  4:  Calling  sequence  for  a  single-point  integration  that  allows  access  to  intermediate  values. 


practical  aspects  of  using  CHEMEQ2. 

Figure  5  illustrates  the  use  of  CHEMEQ2  in  a  reacting  flow  code.  The  effects  of  the  fluid  dynamics  are 
calculated  by  a  separate  algorithm,  and  then  the  conditions  in  each  compuational  cell  are  sent  to  CHEMEQ2 
separately  to  calculate  the  effects  of  the  chemistry.  This  is  the  simplest  implementation  for  reacting  flow, 
and  as  mentioned  earlier  massively  parallel  versions  of  CHEMEQ  have  been  implemented  [14,15].  In  general 
the  optional  calls  to  CHEMSP  will  not  be  needed  since  the  solution  parameters  set  by  CHEMSP  will  be 
the  same  for  all  cells.  If  CHEMSP  must  be  called  repeatedly,  the  logical  check  that  determines  whether 
information  about  the  call  to  CHEMSP  is  printed  can  be  removed. 
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Calculate  initial  conditions, 
control  parameters 


Call  CHEMSP 


advance  fluid  dynamics  At  using 
fluid  dynamics  routines 


optional  CALL  CHEMSP 


Call  CHEMEQ2forcelM 


t  =  t  +  At 


Figure,  5:  Calling  sequence  in  a  reacting-flow  program. 


The  parameter  i-max  gives  the  number  of  cells 


5,4  Practical  Considerations 


To  help  the  beginning  user  of  CHEMEQ2,  some  observations  and  suggestions  are  made  for  optimizing  the 
use  of  the  subroutine.  First,  the  improvement  produced  by  adding  corrector  iterations  is  problem  specific. 
Some  systems,  such  as  cesium  integration  discussed  in  Section  6.1,  converge  in  about  three  iterations.  Other 
systems,  such  as  the  hydrogen- air  mechanism  discussed  in  Section  6.2,  take  much  longer.  Lowering  e  may 
be  more  effective  in  improving  accuracy  than  increasing  Nc  depending  upon  the  problem.  If  a  very  accurate 
result  is  required  and  the  added  computational  cost  can  be  tolerated,  increasing  the  number  of  corrector 
iterations  dramatically  is  very  effective  [19]. 

As  stated  earlier,  a-QSS  is  not  guaranteed  to  be  stable  for  nonlinear  systems.  If  instability  is  seen,  the 
user  can  use  the  convergence-based  stability  check  on  discussed  in  Section  6.2  if  Nc  >  3.  The  lines  which 
implement  this  check  have  a  in  the  first  column.  Many  compilers  allow  the  user  to  include  these  lines 
in  the  compilation  by  specifying  a  compiler  option  such  as  -d^lines.  Without  such  an  option,  these  lines 
are  treated  as  comments  and  not  compiled.  Since  most  reacting-flow  problems  will  not  require  the  stability 
check,  this  implementation  is  most  efficient.  If  the  stability  check  is  needed  on  a  platform  that  does  not 
support  such  a  compiler  option,  then  the  lines  must  be  manually  included. 

The  two  options  for  approximating  a  described  in  Section  4.1  are  included  in  the  code,  but  the  one  labeled 
Fade  (b)  is  recommended.  Although  the  combination  of  Fade  (a)  in  Eq.  (41)  and  the  linear  approximation 
is  Eq.  (39)  are  closer  to  the  exact  curve  for  a  (see  Fig.  2),  the  single  equation  given  in  Eq.  (43)  provides 
comparable  accuracy  in  the  species  concentrations  and  eliminates  an  expensive  logic  check. 

Finally,  users  may  find  the  variable  groupings  in  the  argument  lists  of  CHEMEQ2  and  CHEMSP  incon¬ 
venient.  The  user  may  wish  to  move  arguments  from  one  list  to  the  other,  concentrating  the  parameters 
that  change  regularly  in  the  CHEMEQ2  argument  list,  and  relegating  to  CHEMSP  those  parameters  that 
need  not  be  changed  after  an  initialization.  This  could  save  the  additional  calls  to  CHEMSP  to  reset  a  single 
variable.  For  instance,  if  the  source  term  calculation  depends  upon  the  value  of  the  independent  variable 
then  this  variable  could  be  passed  through  the  CHEMEQ2  argument  list  rather  than  set  through  a  call  to 
CHEMSP. 
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Table  1:  Cesium  Mechanism 


Reaction 

ki 

1) 

0  J  +  Cs+  tlx  CS  +  O2 

5  X  10-s  cmVs 

2) 

Cs+  +  eh  Cs 

1  X  10-12  cm^/s 

3) 

Cs  hCs+  +  e 

3.24  X  10-3  g-i 

4) 

0^h02+€ 

4  X  10-1  s-i 

5) 

02  +  Cs  +  Mh  CSO2  +  M 

1  X  10-31  cm®/s 

6) 

02  +  6  +  02^0^  +O2 

1.24  X  10-30  cmVs 

7) 

02  +  e  +  N2^02  +  N2 

1  X  10-31  cm®/s 

6  Numerical  Results 

Two  examples  are  described  here  in  detail.  The  first  is  a  system  of  equations  involving  cesium  and  cesium 
ions  that  was  originally  suggested  by  D.  Edelson  of  Bell  Laboratories.  This  test  was  used  to  compare  the 
original  CHEMEQ  subroutine  to  other  stiff  solvers,  including  those  of  Gear  and  Kregel,  as  shown  in  [2].  The 
second  set  of  tests  involves  a  hydrogen-oxygen  combustion  mechanism  and  focuses  on  the  effect  of  corrector 

iteration  on  the  timing  and  accuracy  of  a-QSS.  Two  reacting-flow  applications  are  then  discussed  briefly  in 
Section  6.3. 

6.1  Cesium  Tests 

The  cesium  mechanism,  shown  in  Table  1,  involves  seven  species  and  seven  one-way  reactions.  The  rate 
constants  1,  are  fixed  at  the  values  shown.  The  inert  collision  partner,  M  in  reaction  5,  may  be  Cs,  GsOz, 
O2,  or  N2,  so  the  concentration  of  M  used  to  calculate  the  reaction  rate  is  the  sum  of  the  concentrations  of 
these  four  species.  The  initial  conditions  and  the  solution  values  at  1000  seconds  used  for  the  accuracy  study 
are  included  in  the  Table  2  [2].  These  solution  values,  which  we  call  the  “accepted  values”  in  the  following 
error  analysis,  are  the  common  result  of  running  LSODE  and  CHEMEQ2  at  excessively  high  accuracies. 

The  species  number  densities,  shown  in  Fig.  6,  were  generated  using  CHEMEQ2.  The  figure  shows  a 
fast  initial  transient,  which  is  followed  by  a  slow  evolution  toward  equilibrium.  A  logarithmic  scale  in  time 
is  required  to  show  this  evolution.  Equilibrium  is  not  reached  by  1000  s,  so  comparing  the  solution  at  this 
time  to  an  accepted  solution  provides  a  suitable  check  of  a  kinetic  integrator. 
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Number  Density  (cm  *’) 


Table  2:  Initial  and  t  —  1000  sec  species  concentrations  for  the  the  cesium  mechanism  test  problem 


Species 

yi(0s)(cm“^) 

t/i(1000s)[cm 

e 

OJ 

Cs+ 

Cs 

CsOa 

N2 

O2 

1  X  10^ 

5.2  X  10^ 

6.2  X  102 

1  X  10^2 

0 

1.4  X  10‘^ 

3.6  X  lO^'* 

4.9657897283  x  lO'* 
2.5913949444  X  10'* 
7.5571846728  x  10^ 
1.5319405460  x  10^ 
1.000  X  10^2 

1.400  X  lO^® 

3.590  X  10^^ 

Figure  6:  Species  number  densities  as  a  function  of  time  for  the  cesium  mechanism  test  problem 


integration  using  CHEMEQ  and  CHEMEQ2  for  £  =  0.01  and 
1  y  ra  lo,  R,  defined  by  Eq.  (55),  for  O2  ■  R  is  calculated  with  At  used  by  CHEMEQ2. 


Timestep  histories  for  CHEMEQ2  and  for  CHEMEQ  are  shown  in  Fig.  7.  As  mentioned  earlier,  the 

asymptotic  update  used  by  CHEMEQ  is  unstable  under  some  circumstances  [19].  The  linear  stability  analysis 
of  CHEMEQ  led  to  a  parameter  R,  defined  by 


which  we  call  the  stability  ratio.  The  average  timescale,  r,  is  given  by 


(55) 


1^1/1  1  \ 

r-2Vr“  +  ?Fj  (56) 

for  the  initial  value  r°  and  the  predicted  value  r'*.  Stability  requires  R<2  for  any  species  integrated  with 
the  asymptotic  method  [19].  If  the  timescale  is  constant  or  decreasing,  this  stability  constraint  is  satisfied. 
If  the  timescale  is  increasing.  At  may  be  large  enough  to  make  the  method  unstable. 
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During  the  first  200  seconds  of  the  simulation  (i.e.,  the  span  of  time  shown  in  Fig.  7),  CHEMEQ  treats 
only  O2  with  the  asymptotic  update,  so  R  for  this  species  is  included  in  the  figure.  The  values  of  At  are 
read  from  the  vertical  axis  to  the  left  of  the  figure,  and  the  values  of  R  correspond  to  the  axis  on  the  right. 
The  stability  limit  of  /i  =  2  is  marked  by  a  dashed  horizontal  line.  We  see  that  R  becomes  larger  than  2  after 
approximately  10  s,  and  CHEMEQ  starts  producing  oscillations  in  At  approximately  10  s  after  that.  These 
oscillations  cease  after  R  returns  to  values  lower  than  2.  CHEMEQ2  does  not  produce  these  oscillations, 
although  the  accuracy-based  timestep  constraint  lowers  the  timestep  in  this  region. 

A  series  of  studies  evaluated  the  accuracy  of  CHEMEQ2  compared  to  CHEMEQ.  These  solved  the  Cs 
test  problem  given  above  and  used  the  reference  solution  at  1000  s  as  a  benchmark.  The  tests  varied  the 
value  of  €  from  10“^  to  10”^.  Additional  tests  fixed  e  and  varied  Nc  from  one  to  ten.  Figure  8  summarizes 
the  results  of  the  tests  by  showing  the  rms  error  as  a  function  of  CPU  time,  which  was  scaled  by  the  smallest 
increment  the  timing  routine  could  resolve.  The  CHEMEQ2  results  are  shown  as  a  series  of  overlapping 
profiles  of  the  shape  shown  in  the  schematic  in  Figure  9.  Each  profile  is  for  a  fixed  value  of  e,  but  the  points 
on  it  correspond  to  different  values  of  Ac. 

The  error  computed  for  each  computation  (fixed  c  and  Nc)  is  based  on  the  the  accepted  values  at  1000  s. 
The  relative  error  e,-  for  each  species  i 

accepted  ~  2/,*  calculated 

a  = - ^ - .  (57) 

accepted 

A  root-mean-square  error  for  the  six  reacting  species  (excluding  the  inert  N2)  is: 

(58) 

There  is  only  a  single  curve  for  CHEMEQ  in  Figure  5,  Each  point  on  this  curve  corresponds  to  a  different 
value  of  e.  The  hybrid  method,  as  implemented  in  CHEMEQ  and  used  in  this  problem,  becomes  unstable 
and  the  solutions  are  corrupted  if  multiple  corrector  iterations  are  used,  Lorenzini  and  Passoni,  however, 
were  able  to  use  multiple  corrector  iterations  successfully  in  other  implementations  of  the  hybrid  method  for 
other  problems  [29].  CHEMEQ2  does  not  have  this  instability  problem. 

For  a  single-iteration  and  large  enough  c,  the  CHEMEQ2  results  lie  roughly  along  the  CHEMEQ  curve. 
In  this  case,  the  CHEMEQ2  simulation  takes  less  time,  but  gives  a  slightly  less  accurate  solution.  As  e  is 
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Figure  8:  Root  mean  square  error  at  1000  s  versus  scaled 
CHEMEQ2  for  a  range  of  e  and  AT,.  A  schematic  of  each 
numbers  next  to  each  CHEMEQ2  curve  indicate  the  value 
correspond  to  e  =  lO"!  (highest  point  on  the  curve)  to  lO"® 


CPU  time  to  reach  1000  s  for  CHEMEQ  and 
CHEMEQ2  curve  is  shown  in  Fig.  9,  and  the 
^  for  those  results.  The  CHEMEQ  results 
(lowest  point). 


28 


log^Q(scaled  CPU  time) 


Figure  9:  Schematic  of  the  types  of  profiles  for  fixed  £  in  Fig.  8,  The  numbers  next  to  each  symbol  give  the 
corresponding  value  of  Nc^  As  described  in  the  text,  the  solutions  for  the  cesium  test  problem  converge  after 
about  three  iterations. 
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decreased,  the  CHEMEQ  results  do  not  give  the  same  increase  in  accuracy  for  the  increased  computational 
costs. 

The  curves  shown  in  Figure  5  can  be  explained  by  comparing  the  CHEMEQ  and  the  CHEMEQ2  al¬ 
gorithms.  The  CHEMEQ  stiff  predictor  is  identical  to  the  CHEMEQ2  predictor  in  the  limit  as  or®  1, 
which  corresponds  to  p^At  oo.  The  CHEMEQ  stiff  corrector,  however,  uses  different  average  values  for 
q  and  p  than  the  CHEMEQ2  corrector,  and  also  effectively  uses  the  pAt  ^  0  limit  value  of  a  =  0.5.  This 
inconsistency  in  the  effective  a  between  CHEMEQ’s  stiff  predictor  and  corrector  limits  the  growth  of  At  for 
CHEMEQ.  Therefore,  CHEMEQ  takes  a  smaller  timestep  than  CHEMEQ2  for  the  same  e,  and,  for  moderate 
accuracy,  this  inconsistency  in  a  does  not  affect  the  accuracy  of  the  solution.  The  best  accuracy  achievable 

by  CHEMEQ  does  suffer  from  this  inconsistency,  however,  so  as  e  becomes  smaller,  CHEMEQ2  gives  more 
accurate  answers  than  CHEMEQ. 

The  CHEMEQ2  curves  for  a  fixed  e  show  dramatic  increases  in  accuracy  after  just  a  few  iterations.  After 
about  three  iterations,  the  curves  for  a  given  e  flatten,  which  indicates  that  the  method  has  converged  to  a 
final  corrector  value,  and  additional  iterations  do  not  improve  the  accuracy.  The  computational  expense  in 
adding  iterations  is  less  than  that  in  reducing  e  for  similar  improvements  in  accuracy.  As  e  is  lowered, 
accuracy  improves  because  the  timestep  is  decreased.  As  the  number  of  iterations  increases,  accuracy 
improves  because  the  corrector  is  able  to  refine  the  linear  approximation  for  p  and  q  used  to  calculate  q 
and  p  for  the  corrector  equation,  Eq.  (36).  Not  all  systems  will  converge  for  such  low  values  of  A,,  but,  in 
general,  iterating  the  corrector  even  one  or  two  times  improves  the  accuracy. 

For  the  CHEMEQ2  curve  for  £  =  0.1,  the  simulations  took  so  little  time  that  the  precision  of  the  timing 
routines  was  not  sufficient  to  measure  differences  in  timing  between  these  runs.  In  addition,  the  calculations 

were  performed  on  a  computer  that  allows  access  to  multiple  users.  These  affects  contribute  to  the  error  and 
uncertainty  in  the  low-resolution  data. 

6.2  Hydrogen- Air  Tests 

The  Ao-air  combustion  mechanism  used  consists  of  twenty-five  reversible  reactions  involving  nine  species 
(including  inert  A,)  [34].  This  mechanism  is  clo.sely  related  to  that  used  by  GRIMech.  The  reaction  rales 
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Figure  10:  Solution  for  the  single-point  hydrogen-air  integration. 

are  calculated  using  the  modified  Arrhenius  form 

kr  =  AT^cxp{-C/T)  (59) 

where  T  is  the  temperature.  The  rate  kr  is  either  a  forward  or  backward  rate.  The  parameters  A,  B,  and  C 
for  both  the  forward  and  backward  rates  for  each  reaction  are  given  in  reference  [19].  Initially  the  mixture 
is  at  1000  /i,  a  pressure  of  1  atmosphere,  and  in  the  ratio  2:1:3.76  for  These  conditions  lead 

to  initial  number  densities  on  the  order  of  lO^^cm”^  for  these  three  species.  A  minimum  number  density 
of  10"^^cm”^  wcis  imposed  on  the  other  species  to  prevent  numerical  difficulties.  Nitrogen  is  inert  for  the 
mechanism,  and  thus  acts  as  a  diluent. 

Selected  species’  number  densities  for  this  problem  are  presented  as  a  function  of  time  in  Fig.  10.  The 
figure  shows  that  after  an  induction  time  of  about  3.4  x  10”"^  s,  Hn  and  On  are  converted  to  HnO  in  a 
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relatively  short  time  period.  During  this  induction  time  radicals  are  formed  that  eventually  initiate  the 
rapid  conversion  of  and  Oj.  Here,  we  focus  on  the  H  number  density  profile,  which  has  a  peak  in  the 
reaction  zone  that  is  difficult  to  predict  accurately.  A  series  of  calculations  examined  the  effect  of  e  and  N, 
on  the  location  and  the  value  of  this  peak.  The  errors  in  these  parameters  are  calculated  as 

i  ~  ^P.referencell 

p  ; -  (60) 

P.reterence  ^ 


(ni/)p  error  =  — reference II 


, reference 


for  peak  number  density  value  (nff)p  at  time  tp,  compared  to  reference  values.  The  reference  values  were 
obtained  by  integrating  the  equations  with  CHEMEQ2  for  increasing  AT,  and  decreasing  s  values,  until 
the  solution  ceased  changing.  The  solution  was  then  verified  by  comparison  with  a  solution  obtained  by 
a  simple  modified  Euler  method  using  an  exceptionally  low  error  tolerance.  Table  3  lists  these  errors  and 
the  CPU  time  required  to  reach  5  x  lO"'*  seconds  for  a  variety  of  s  and  N.  values.  These  calculations  were 
performed  on  a  DEC  Alpha  workstation,  and  the  integration  time  is  scaled  relative  to  the  s  =  =  1 

simulation.  These  calculations  did  not  assume  that  the  thermodynamic  state  or  the  rate  constants  remained 
fixed  during  a  chemical  timestep.  The  temperature  was  recalculated  for  each  corrector  iteration  based  on  the 
species  number  densities  calculated  from  the  previous  iteration.  The  rate  constants  were  then  recalculated 
with  this  updated  temperature.  Repeating  this  calculation  for  every  iteration  is  very  expensive,  and  the 

performace  of  the  method  will  improve  substantially  if  the  rate  constants  are  calculated  once  and  fixed  for 
the  chemical  timestep. 


Figure  11  shows  results  of  integrations  for  e  =  10"^  and  N,  =  1,  5,  and  10.  This  should  be  contrasted 
to  the  cesium  calculations  of  the  previous  section  that  converged  by  N,  =  5.  In  this  case,  the  profiles 
are  converging  to  the  reference  solution,  but  they  have  not  completely  converged  by  N,  =  10.  Note  that 
CHEMEQ  results  are  essentially  equivalent  to  the  N,  =  1  case.  Table  3  suggests  that  reducing  e  may  be  a 
more  efficient  way  to  improve  the  accuracy  of  the  solution  than  increasing  the  iteration  count.  The  errors  are 


not  of  the  same  order  as  e,  however,  and  reducing  s  by  an  order  of  magnitude  does  not  result  in  a  comparable 
reduction  in  the  error.  The  errors  in  the  time-to-peak  and  the  peak  value  are  not  even  comparable,  with  the 
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Figure  11:  Hydrogen  number  density  for  Nc  =  1,  5,  and  10,  and  e  =  10  The  dark,  solid  line  is  the 
reference  solution,  and  the  numbers  next  to  the  remaining  curves  indicate  the  value  of  Nc  for  each  profile. 
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Table  3:  Results  obtained  by  varying  e  and  TV,  for  the  hydrogen-air  reaction  integration. 


3(a):  tp  error 


€ 

Nc  =  l 

5 

10 

10-3 

6.66  X  10“’^ 

2.97  X  10-3 

”1.84  X  10-3 

10"^ 

2.79  X  10-3 

1.10  X  10-3 

6.29  X  10-3 

10-3 

1.06  X  10-3 

3.67  X  10-3 

1.97  X  10-3 

3(b): 

error 

£ 

Nc  =  l 

5 

10 

10-3 

0.392 

0.166 

9.40  X  10-3 

10"^ 

0.146 

4.56  X  10-3 

2.35  X  10-3 

10-3 

4.48  X  10-3 

1.27  X  10-3 

6.49  X  10-3 

3 

- - - - — - - - - 1 

[c):  Scaled  CPU  times  to  5  x  lO""*  s. 

e 

Nc  =  l 

5 

"  10 

0 

1 

o 

r-^ 

1.00 

2.92 

5.33 

10-^ 

3.19 

9.92 

18.3 

o 

1 

Cn 

11.8 

36.7 

67.5 

Table  4:  Errors  in  tp  and  (n„)p  for  c  =  IQ-^  and 
for  each  simulation  to  reach  <  =  5  x  10"^  seconds. 


1,  5,  10,  and  1000,  and  the  scaled  CPU  time  required 


iterations 

tp  error 

Ph  error 

CPU  time 

1 

6.66  X  10-3 

0.392 

1.00 

5 

2.97  X  10-3 

0.166 

2.92 

10 

1.84  X  10-3 

9.40  X  10-3 

5.33 

1000 

2.71  X  10-3 

1.77  X  lO-'* 

489 

Mp  much  more  prone  to  error  than  tp.  This  peak  is  very  difficult  for  a  low-order  method  to  calculate.  A 

higher-order  method  that  employs  information  from  several  timesteps  would  provide  a  much  better  result 
for  this  problem. 

The  question  remains  as  to  how  accurate  the  integration  can  become  if  the  number  of  iterations  is 
increased  dramatically.  The  results  for  e  =  10-3  from  Table  3  are  repeated  in  Table  4,  and  additional  results 
obtained  using  1000  corrector  iterations  are  also  included.  For  the  N,  =  1000,  the  error  in  the  peak  value 
is  an  order  of  magnitude  less  than  e,  and  the  time-to-peak  error  is  three  orders  of  magnitude  lower  than  s. 
This  suggests  that  the  corrector  equation  Eq.  (36)  provides  an  accurate  representation  once  it  is  sufficiently 
converged.  Again,  the  thermodynamic  state  and  the  rate  constants  were  recalculated  for  every  corrector 
Iteration,  and  the  high  CPU  times  are  due  in  part  to  this  largely  avoidable  c.xpense. 
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Figure  12:  Hydrogen  number  density  as  the  integration  approaches  equilibrium,  e  —  10“'^,  Nc  ~  10.  The 
dashed  line  is  the  standard  CHEMEQ2  result.  The  profile  given  by  open  circles  includes  the  stability 
constraint  on  A<  (see  Eq.  (52)), 

An  instance  which  requires  the  stability  check  on  Ai  described  in  Section  4.2  is  illustrated  in  Fig.  12.  This 
figure  shows  the  H  profile  as  the  system  reaches  equilibrium.  These  results  indicate  that  the  accuracy-based 
timestep  can  be  too  large  for  the  corrector  iteration  to  remain  stable  despite  the  fact  that  the  stability  analysis 
indicated  that  a-QSS  is  A-stable  for  linear  problems.  Note  that  the  scale  in  Fig.  12  is  exaggerated;  the  range 
covered  by  the  number  density  axis  spans  approximately  1%  of  the  equilibrium  value.  This  instability  is  not 
a  problem  in  reacting-flow  applications,  as  the  frequent  restarting  at  new  global  timesteps  limits  how  large 
the  timestep  becomes.  In  this  single-point  integration,  however,  the  instability  is  seen.  The  oscillations  in 
the  number  density  disappear  when  the  stability  constraint  given  in  Eq.  (52)  is  required,  and  the  predicted 
equilibrium  value  agrees  well  with  the  reference  solution. 
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6.3  Reacting-FIow  Solutions 


Two  reacling-aow  casas  will  be  briefly  discussed  here.  These  results  are  pro.isional,  as  no  rigorous,  systematic 
studies  ha.e  been  performed.  A  thorough  comparison  between  integration  methods  would  include  the  effects 
of  .mplementation  choices,  accuracy  requirements,  and  stiffness.  The  stiffness  issues  ate  not  limited  to 
the  chemical  mechanism  itself  but  also  include  coupling  of  the  chemical  timescales  and  the  fluid  dynamics 
timescales  (i.e.,  how  much  the  integrator  subdivides  tbe  global  timestep  in  order  to  perform  the  chemistry 

integration).  Such  a  study  is  planned  for  the  future.  However,  from  out  experiences,  we  expect  the  results 
described  below  to  be  typical. 

Uphoff  et  al.  [35]  studied  two-dimensional  detonation  formation  using  an  mechanism  with  18  re¬ 

actions  and  8  species.  They  compared  process-split  reacting-flow  calculations  using  CHEMEQ  and  METANl 
[36]  as  the  chemistry  integrator.  METANl  is  a  general  stiff  solver  which  employs  a  semi-implicit  mid-point 
rule  and  extrapolation  to  a  “zero  stepsize”  solution  [37-39].  For  this  specific  set  of  calculations,  CHEMEQ 
performed  the  required  chemical  integrations  in  approximately  one-sixth  of  the  time  required  by  METANl. 
Documentation  of  accuracy  parameters  used  and  solution  options  chosen  for  the  calculations  is  not  available. 

An  additional  calculation  was  performed  in  order  to  compare  the  efficiency  of  o-QSS  to  a  Gear  method.  A 
one  dimensional  hydrogen-air  premixed  flame  was  simulated  using  a  process-split  method  [33]  which  employed 
FCT  for  integrating  the  fluid  convection  [40].  The  chemistry  integration  was  performed  using  CHEMEQ2, 
and  also  using  DEBDF,  which  employs  a  variable-order  Gear  method  as  implemented  in  LSODE.  DEBDF  is 
part  of  SLATEC,  a  library  of  computational  subroutines  available  on  Silicon  Graphics  and  Cray  computers 
[41].  CHEMEQ  performed  the  required  calculations  in  approximately  one-sixth  the  time  required  by  DEBDF, 
which  is  coincidentally  the  same  factor  seen  in  the  detonation  comparison  versus  METANl.  No  extensive 
accuracy  studies  have  been  performed  to  ensure  that  the  comparison  was  fair.  For  example,  the  accuracy 

parameters  for  CHEMEQ2  and  DEBDF  were  simply  set  to  the  same  value,  even  though  the  two  codes  do 
not  use  these  parameters  in  exactly  the  same  way. 
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7  Summary 


CHEMEQ2  is  a  general  purpose  integrator  for  a  specific  type  of  equations,  namely  those  that  are  reasonably 
represented  by  the  form  in  Eq.  (1).  CHEMEQ2  employs  a  very  low  overhead,  moderately  accurate,  low-order 
technique.  To  obtain  results  for  most  physical  models  with  an  acceptable  degree  of  accuracy,  CHEMEQ2  can 
be  extremely  efficient.  In  many  areas  where  problems  are  so  computationally  expensive  they  seem  impossible 
to  do  by  other  methods,  CHEMEQ2  gives  accurate  results  in  a  reasonable  amount  of  time.  CHEMEQ2  can 
also  be  employed  in  the  development  of  chemical  or  mathematical  models  when  efficiency  is  important,  but 
obtaining  very  precise  answers  may  require  extensive  computational  expense.  CHEMEQ2  is  optimized  to 
provide  three  or  four  significant  digits  accurately,  not  eight,  but  this  high  level  of  accuracy  can  be  reached 
with  an  appropriate  timestep  criterion  and  enough  corrector  iterations. 

CHEMEQ2’s  forte  lies  in  the  solution  of  the  stiff  ordinary  differential  equations  associated  with  chemically 
reactive  flow  problems.  Here  the  reaction  sources  are  split  off  from  the  hydrodynamic  part  of  the  equations 
and  solved  separately  for  each  hydrodynamic  timestep  and  at  each  grid  point.  The  moderate  accuracy  of 
the  methods  used  to  solve  the  hydrodynamic  equations  suggest  that  the  application  of  a  more  sophisticated 
technique,  rather  than  a  low-order,  low  overhead  method  like  CHEMEQ2,  would  waste  valuable  computer 
time  and  could  possibly  render  the  problem  impossible. 

A  potential  user  must  be  aware  that  CHEMEQ2  is  not  completely  user-proof  or  problem-independent 
and  cannot  always  be  used  as  a  black  box.  The  method  is  not  identically  conservative,  and  the  minimum 
values  should  be  chosen  with  some  thought  since  they  can  become  sources  of  spurious  errors  if  not  chosen 
small  enough  initially.  Although  CHEMEQ2  overcomes  some  stability  problems  in  the  original  algorithm, 
it  may  still  require  the  use  of  the  stability  check  described  in  Section  4.2.  The  user  is  referred  again  to 
Section  5.4  for  practical  information  regarding  the  use  of  CHEMEQ2. 

Since  CHEMEQ2  uses  a  convergence-dependent  algorithm  and  an  adaptive  timestep,  the  overall  timing 
will  be  strictly  problem-dependent.  One  factor  will  be  the  coupling  between  the  relaxation  times  of  the 
equations.  The  most  expensive  operation  in  the  algorithm  is  the  derivative  function  evaluations,  of  which 
there  is  one  required  in  the  predictor  step,  and  one  for  each  corrector  iteration.  If  CriEMEQ2  is  applied  as 


designed,  the  subroutine  can  solve  large  systems  of  stiff  ordinary  differential  equations  very  efficiently. 
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Appendix  A  -  CHEMEQ2  Variable  List 

Table  9:  Variable  listing  and  descriptions. 


FORTRAN 

variable 


Type/Scope  Same  As 


Description 


alpha 

R/L 

or,  Eq.  (7) 

ascr 

R/L 

— 

d(i) 

R/A 

pm,  Eq.  (1) 

dt 

R/L 

At 

dto 

R/L 

— 

dtc 

R/L 

mhi 

dtg 

R/A 

Atg 

dtmin 

R/L 

— 

dtmn 

R/A 

— 

epscl 

R/L 

1/epsmin 

eps 

R/L 

<r,  Eq.  (48) 

epsmax 

R/A 

c  from  Eq.  (47) 

epsmin 

R/A 

£  from  Eq.  (46) 

epsmn 

R/A 

_ 

epsmx 

R/A 

— 

gcount 

I/L 

— 

gsub 

E/A 

_ 

i 

I/L 

i 

iter 

I/L 

— 

itermax 

I/A 

— 

itermx 

R/A 

— 

lo 

I/L 

— 

n 

I/A 

n 

nd 

I/L 

— 

ns 

I/A 

— 

prt 

R/A 

— 

q(i) 

R/A 

9t,  Eq.  (1) 

qs(i) 

R/A 

9?,  Eqs.  (35) 
(38) 

qt 

R/A 

9£,  Eq.  (38) 

rcount 

I/L 

— 

Type:  R  =  Real,  I  = 


solution  parameter  used  in  update 
scratch  (temporary)  variable 
loss  rate 

chemical  timestep  used  by  the  integrator 
stores  timestep;  used  to  scale  rtaus  when  timestep 
is  reduced 

diagnostic  value  printed  when  At  <  dtmin 
global  timestep;  range  of  integration 
minimum  timestep  allowed 
sets  dtmin  via  CHEMSP 

intermediate  variable  used  to  avoid  repeated 
divisions 

maximum  correction  term,  finally  scaled  by 
1/epsmin 

repeat  timestep  if  correction  is  greater  than 
epsmax*epsmin*y(i)  for  any  i 

accuracy  parameter  for  determining  the  next 
timestep 

used  to  set  epsmin  via  CHEMSP 
used  to  set  epsmax  via  CHEMSP 

counter  for  calls  to  gsub  since  the  last  call  to 
CHEMCT 

source  term  subroutine;  supplies  d(i)  and  q(i) 
index 

counter  for  corrector  iterations 
number  of  corrector  iterations  to  perform 
used  to  set  itermax  via  CHEMSP 
unit  number  for  output 
number  of  equations  integrated 

dimension  of  species  arrays;  maximum  number  of 
species 

number  of  entries  in  ymin  reset  via  CHEMSP  call 
nonzero  value  supresses  output  from  CHEMSP 
production  rate 
initial  production  rate 

a- weigh  ted  average  of  q 

counter  for  steps  redone  since  the  last  call  to 
_ _ CHEMCT 

li^eger.  E  =  External;  Scope:  L  ^  Strictly  Local.  A  =  Passed  as  Argument 
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Table  9  Continued 


FORTRAN 

variable 

Type/Scope 

Same  As 

Description 

rswitch 

R/L 

5.9659 

value  of  A^/r  used  to  switch  between  Eqs. (39)  and 
(41)  when  Fade  (a)  is  used 

rtau(i) 

R/L 

At/n 

ratio  of  timestep  to  timescale 

rtaub 

R/L 

piAt  = 

At  In 

times  average  p  from  Eq.  (37) 

rtaui 

R/L 

Ai/Ti 

holds  rtau(i)  to  avoid  multiple  array  references 

rtaus(i) 

R/L 

At/rf 

ratio  of  timestep  to  initial  timescale  for  current 
timestep 

rteps 

R/L 

estimate  for  y/a  in  Eq.  (50) 

scrl 

R/L 

— 

scratch  (temporary)  variable 

scr2 

R/L 

— 

scratch  (temporary)  variable 

scraxray 

R/L 

— 

scratch  (temporary)  variable  array 

scrtch 

R/L 

— 

scratch  (temporary)  variable 

sqreps 

R/A 

by/e 

parameter  used  to  calculate  initial  timestep 

stab 

R/L 

— 

see  Eqs.  52)  and  (53) 

tfd 

R/L 

round-off  parameter  used  to  determine  when  inte¬ 
gration  is  complete 

tgcnt 

I/L 

total  number  of  calls  to  gsub  for  all  global 
timesteps 

tmk 

R/A 

— 

call  identifier  for  CHEMCT 

tn 

R/A 

t-t^ 

current  value  of  the  independent  variable  relative 
to  the  start  of  the  global  timestep 

tnot 

R/A 

— 

used  to  set  tstart  via  CHEMSP 

trcnt 

I/L 

— 

total  number  of  steps  redone  for  all  global 
timesteps 

ts 

R/A 

independent  variable  at  the  start  of  the  chemical 
timestep 

t start 

R/A 

independent  variable  at  the  start  of  the  global 
timestep 

y(i) 

R/A 

Vi 

species  concentrations  array 

yO(i) 

R/A 

from  Eq.  (2) 

initial  concentrations  for  the  global  timestep 
passed  to  CHEMEQ 

yl{i) 

R/A 

predicted  value  from  Eq.  (35) 

yml(i) 

R/L 

previous  corrector  iterate;  see  Eq.  (51) 

ym2(i) 

R/L 

previous  corrector  iterate;  see  Eq.  (51) 

yinin(i) 

R/L 

— 

minimum  concentration  allowed  for  species  i 

yinin(i) 

R/A 

— 

set  ymin(i)  via  CHEMSP 

ys(i) 

R/L 

Eqs. 

(36) 

from 
(35)  and 

initial  concentrations  for  the  chemical  timestep 

Type:  R  = 

Real,  I  =  Integer,  E  = 

External;  Scope:  L  =  Strictly  Local,  A  =  Passed  as  Argument 
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10 

10.1 

c 

cd*  * 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 


Appendix  B:  Code  Listings 
CHEMEQ2  Code  Listing 

subroutine  chemeq2(dtg,  gsub,  n,  y) 

cheineq2(dtg,  gsub,  n,  y) 


original  chemeq  development : 
originators:  t.r.  young 
vax  version:  t.r.  young 
workstation:  g.  patnaik 

chemeq2  development:  d.r.  mott 


nrl  1982 

nrl  code  4040  may  1983 

berkeley  research  jun  1995 

nrl  code  6404  may  1999 


Description:  Subroutine  chemeq2  solves  a  class  of  ^stiff’  DDEs 
associated  with  reactive  flow  problems  that  cannot  be  readily 
solved  by  the  standard  classical  methods.  In  contrast  to  the 
original  chemeq  subroutine,  this  version  uses  the  same 
quasi-steady-state  update  for  every  species  regardless  of  the 
timescale  for  that  species.  An  adaptive  stepsize  is  chosen  to 
give  accurate  results  for  the  fastest  changing  quantity,  aud  a 
stability  check  on  the  timestep  is  also  available  when  the 
corrector  is  iterated. 

NOTE;  The  accuracy-based  timestep  calculation  can  be  augmented 
with  a  stability-based  check  when  at  least  three  corrector 
iterations  are  performed.  To  include  this  check,  ••uncomment" 
the  lines  that  start  with  '•D^^,  or  use  the  compiler  flag  "-d_lines^^ 
if  available  to  compile  the  code  including  these  lines.  If  the 
lines  are  manually  uncommented,  the  continuation  characters 
must  be  placed  in  the  correct  column.  For  most  problems,  the 
stability  check  is  not  needed,  and  eliminating  the  calculations 
and  logic  associated  with  the  check  enheinces  performance. 

The  routine  assumes  that  all  of  the  integrated  quantites  cind  the 
time  step  are  positive. 


argument  list  definition 
real 

gsub  real 

^  integer 


y(n)  real 


(name,  type,  description,  input  vs.  output): 
the  interval  of  integration  or  the  i 
range  of  the  independent  variable, 

0.0  <=  t  <=  dtg.  (global  timestep) 
the  name  of  the  derivitive  function  i 
evaluator  subroutine. 

the  number  of  equations  to  be  i 

integrated,  an  error  exisis  if  n  is 
greater  than  nd  set  by  the  parameter 
statement , 

the  initial  values  at  call  time  i/o 
and  the  final  values  at  return  time. 
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cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd* 

c 


Language  and  limitations:  This  subroutine  is  written  in  standard 
FORTRAN  77.  For  high  accuracy,  this  routine  should  be  compiled 
using  whatever  “double  precision'*  flag  is  appropriate  for  the 
platform  being  used  (such  as  “f77  -r8  .  .  .  .“) 

Entry  points:  Four  entry  points  are  provided  for  flexibility  and 
optimum  control.  This  structure  was  maintained  from  the  original 
chemeq  subroutine  to  ensure  compati ability  with  previous 
applications  that  use  chemeq. 

chemeq2:  advamces  the  equations  the  given  increment  'dtg’. 

chemct:  informative,  prints  the  values  of  the  indicative 
counters  listed  below; 

1,  the  number  of  derivative  function  evaluations. 

2.  the  number  times  the  integration  step  was  restarted 
due  to  nonconvergence  of  the  predictor-corrector 
scheme . 

chemsp:  provides  the  user  with  the  option  to  reset  the  most 
important  control  parameters. 

chempr:  informative,  prints  out  internal  vaxiables  for  diagnostic 
purposes. 

subroutines  referenced: 

gsub;  whose  actual  name  and  definition  are  supplied  by  the  user 
is  called  to  obtain  the  derivitive  functions. 

call  gsub(y,  q,  d,  t) 
argument  list  to  gsub; 

y(n)  real  current  values  of  the  dependent  i 

variable , 

q(n)  real  calculated  formation  rates.  o 

d(n)  real  calculated  loss  rates.  o 

t  real  current  value  of  the  independent  i 

variable . 

chemer:  Called  whenever  an  error  is  detected.  Currently  the 

only  error  recognized  is  a  time  step  that  is  too  small. 

call  chemer (y,  n) 

argument  list  to  chemer;  (same  definition  as  “chemeq2"). 

implicit  none 
integer  nd 

parameter  (nd  =  10) 


c 
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external  gsub 

integer  n,  ns,  lo,  i 

integer  itermax,  iter,  itermx 

the  following  are  counters  (this  call  &  total)  for  gsub  calls 
and  timestep  repeats 


integer  gcount,  rcount,  tgcnt,  trcnt 
c 


c 

c 

D 

c 


real 

real 

real 

real 

real 

real 

real 

real 

real 

real 

real 

real 


ts,  tn,  tfd,  tmk 
y(n) 

ymin(nd),  yinn(nd) 

q(nd),  d(nd),  rtaus(nd),  yl(nd) 

ys(nd),  yO(nd),  rtau(nd) 

alpha,  qs(nd) 

scrl,  scr2,  scrarray(nd) 

epscl,  dtg,  dtmin,  sqreps,  t start,  dt,  dto 

epsmax,  epsrain,  rswitch 

epsmx,  epsmn,  dtmn,  tnot,  prt 

scrtch,  ascr,  eps 

rtaui,  rtaub,  qt,  pb,  dtc,  rteps 


yml,  ym2,  and  stab  are  used  only  for  the  stability  check  on  dt 
real  yml(nd),  yin2(nd),  stab 


data 

data 

data 

data 

data 

data 


gcount,  rcount,  tgcnt,  trcnt /4*0/ 
itermax/l/,  epscl/100.0/ 
tfd/l . 000008/ ,  dtmin/ 1 . Oe-lS/ ,  sqreps/0.50/ 
t St art,  dt/2*0.0/,  tn/0,0e+00/,  q/nd*0.0/ 
epsmax/10 .0/,  lo/16/,  epsmin/1 .Oe-02/,  d/nd*0.0/ 
rswitch/  5.965900  / 


c  rswitch  for  4-4  pade:  5.9659 
c 


cd 


1002 


check  input  parameters , 
if(n  .gt.  nd)  then 

writedo,  1002)  n,  nd 
format (5(/) , ’from  -chemeq2-  : 
'  large’/’  requested  (’,i5,’), 
stop 
end  if 


no.  of  eq.s  requested  is  too’, 
max,  allowed  (’,i5,’)’) 


c 

c  initialize  the  control  parconeters. 
110  tn  =  O.Oe+00 


store  and  limit  to  ’ymin’  the  initial  values 
do  i  =  1,  n 
q(i)  =  0.0 
d(i)  =  0,0 
y0(i)  =  y(i) 

y(i)  =  max(y(i),  ymin(i)) 
end  do 
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c  evaluate  the  derivitives  of  the  initial  values. 

call  gsub(y,  q,  d,  tn  +  tstart) 
gcount  =  gcount  +  1 
c 

c  estimate  the  initial  stepsize. 
c 

c  strongly  increasing  functions (q  »>  d  assumed  here)  use  a  step- 
c  size  estimate  proportional  to  the  step  needed  for  the  function  to 

c  reach  equilibrium  where  as  functions  decreasing  or  in  equilibrium 

c  use  a  stepsize  estimate  directly  proportional  to  the  character- 
c  istic  stepsize  of  the  function,  convergence  of  the  integration 
c  scheme  is  likely  since  the  smallest  estimate  is  chosen  for  the 
c  initial  stepsize. 

scrtch  =  l.Oe-25 
do  i  =  1,  n 

ascr  =  abs(q(i)) 

scr2  =  signd  ./y(i)  ,  .  l*epsmin*ascr  -  d(i)  ) 
scrl  =  scr2  *  d(i) 

scrtch  =  max(scrl,-abs(ascr“d(i))*scr2, scrtch) 
end  do 

dt  =  min(sqreps/scrtch,dtg) 
c 

c  the  starting  values  are  stored. 

100  ts  =  tn 

c 

do  i=l,n 

rtau(i)  =  dt*d(i)/y(i) 
ys(i)  =  y(i) 
qs(i)  =  q(i) 
rtaus(i)  =  rtau(i) 
end  do 
c 
c 

c  find  the  predictor  terms. 

101  continue 
c 

do  i  =  l,n 
c 

c  prediction 

c 

rtaui  =  rtau(i) 
c 

c  note  that  one  of  two  approximations  for  alpha  is  chosen: 

c  1)  Fade  b  for  all  rtaui  (see  supporting  memo  report) 

c  or 

c  2)  Fade  a  for  rtaui<=rswitch , 

c  linear  approximation  for  rtaui  >  rswitch 

c  (again,  see  supporting  NRL  memo  report  (Mott  et  al.,  2000)  ) 

c 

c  Option  1):  Fade  b 


^9 


alpha  -  (180.+rtaui*(60.+rtaui+(H.+rtaui))) 

&  /  (360.  +  rtaui*(60.  +  rtaui*(12.  +  rtaui))) 

Option  2);  Fade  a  or  linear 

if (rtaui. le.rswitch)  then 

alpha  =  (840.+rtaui+(140.+rtaui*(20.+rtaui))) 

*  /  (1680.  +  40.  ♦  rtaui*rtaui) 

else 

alpha  =  l.-l. /rtaui 
end  if 

scrarray(i)  =  (q(i)-d(i))/(l.o  +  alpha*rtaui) 
end  do 

iter  =  1 

do  while(iter .le . itermax) 


limit  decreasing  functions  to  their  minimum  values 
do  i  =  l,n 

3rm2(i)  =  yml(i) 
yml(i)  =  y(i) 

y(i)  =  max(ys(i)  +  dt*scrarray(i) ,  yinin(i)) 
end  do 

if (iter , eq, 1)  then 

the  first  corrector  step  advances  the  time  (tentatively)  and 

saves  the  initial  predictor  value  as  yl  for  the  timestep  check  later 
tn  =  ts  +  dt 
do  i=l,n 

yi(i)  =  y(i) 

end  do 
end  if 


evaluate  the  derivitives  for  the  corrector. 

call  gsub(y,  q,  d.  tn  +  tstart) 
gcount  =  gcount  +  1 
eps  =  l.Oe-10 

do  i  =  l,n 

rtaub  =  .5+(rtaus(i)+dt+d(i)/y(i)) 

Same  options  for  calculating  alpha  as  in  predictor: 

Option  1):  Fade  b 

alpha  -  (180 .+rtaub*(60.+rtaub+(ll .+rtaub) ) ) 

^  /  (360.  +  rtaub* (60.  +  rtaub* (12.  +  rtaub))) 


c 

c 

c 

c 

c 

c 

c 

c 


Option  2):  Fade  a  or  linear 

if (rtaub.le .rswitch)  then 

alpha  =  (840 .+rtaub*(140 .+rtaub*(20.+rtaub) ) ) 
&  /  (1680.  +  40.*rtaub*rtaub) 

else 

alpha  =  l,-l./rtaub 
end  if 


qt  =  qs(i)*(l.  ~  alpha)  +  q(i)*alpha 
pb  =  rtaub/dt 

scrarray(i)  =  (qt  ~  ys(i)*pb)  /  (1.0  +  alpha+rtaub) 
c 

end  do 
c 

iter  =  iter  +  1 
c 

end  do 
c 

c  calculate  new  f,  check  for  convergence,  and  limit  decreasing 
c  functions,  the  order  of  the  operations  in  this  loop  is  important, 
do  i  =  l,n 

scr2  =  max(ys(i)  +  dt*scrarray(i) ,  0,0) 
scrl  =  abs(scr2  -  yl(i)) 
y(i)  =  max(scr2,  ymin(i)) 
ym2(i)  =  3^l(i) 
yml(i)  =  y(i) 

if ( .25*(ys(i)  +  y(i) ) .gt ,ymin(i) )  then 
scrl  =  scrl/y(i) 
eps  =  max(,5*(scrl+ 

min(abs(q(i)“*d(i) )/(q(i)+d(i)+l .Oe-30) , scrl) ) ,eps) 

end  if 
end  do 

eps  =  eps+epscl 
c 

c  print  out  dianostics  if  stepsize  becomes  too  small. 


if(dt  .le.  dtmin  +  1.0e-16*tn)  then 
write(lo,  1003)  dt ,  tn,  dtmin 
do  i  =  l,n 

dtc  =  epsmin*y(i)/(abs(q(i)-d(i) )  +  l.Oe-30) 
write(lo,  1004)  q(i) ,  d(i),  y(i),  rtau(i),  dtc, 
&  q(i)-d(i) ,ys(i) ,  y0(i),  ymin(i) 

end  do 


1003  format(*l  chemeq  error;  stepsize  too  small  !  !  !*,  /, 

1  ’  dt  =  \  lpelO.3,  »  tn  =  \  d25.15, 

2  ’  dtmin  =  ’,el0,3,  //,  llx,  ’q\  lOx,  M’,  lOx,  ’y\ 

3  8x,  ’rtau’,  8x,  Mtc\  7x,  ^q-d\7x,  'ys\ 
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4  9x.  >yO’,  8x,  ’ymin’) 

1004  format(5x,  lpl2ell.3) 

dt  =  dtg  -  ts 

dt  =  min(dtmin,  abs(dt)) 
c 

c  call  error  diagnostic  routine 

call  cheraer 
c 


c 

c 

c 

c 

D 

D 

D 

D 

D 

D 

D 


end  if 

check  for  convergence. 

The  following  section  is  used  for  the  stability  check 
Stab  =0.01 
if (itennax.ge.3)  then 
do  i=l,n 

stab  =  max(stab.  abs(y(i)-yml(i) )/ 

ft  (abs(ynil(i)-ym2(i))+l.e-20*y(i))) 

end  do 
endif 


if(eps  .le.  epsmax 
D  &  .and. stab. le.l. 

&  )  then 

c 

C  Valid  step.  Return  if  dtg  has  been  reached, 
c 

if (dtg  .le.  tn*tfd)  return 
else 
c 

c  Invalid  step;  reset  tn  to  ts 
c 

tn  =  ts 
end  if 
C 

c  perform  stepsize  modifications, 

c  estimate  sqrt(eps)  by  newton  iteration, 

c 

rteps  =  0.5*(eps  +  1.0) 
rteps  =  0.5* (rteps  +  eps/rteps) 
rteps  =  0.5* (rteps  +  eps/rteps) 
c 

dto  =  dt 

dt  =  min(dt*(1.0/rteps  +  .005),  tfd*(dtg  -  tn) 
D  &  ,dto/(stab+.001) 

&  ) 

c 

c  begin  new  step  if  previous  step  converged. 

if(eps  .gt.  epsmax 
D  &  .or.  stab.  gt.  1 

&  )  then 
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rcount  =  rcount  +  1 


c 

c  After  an  unsuccessful  step  the  initial  timescales  don’t 

c  change,  but  dt  does,  requiring  rtaus  to  be  scaled  by  the 

c  ratio  of  the  new  and  old  timesteps. 

c 

dto  =  dt/dto 
i  =  1 

do  while(i . le.n) 

rtaus (i)  =  rtaus (i) ♦dto 
i  =  i+1 
end  do 
c 

c  Unsuccessful  steps  return  to  line  101  so  that  the  initial 

c  source  terms  do  not  get  recalculated, 

c 

goto  101 
end  if 
c 

c  Successful  step;  get  the  source  terms  for  the  next  step 

c  and  continue  back  at  line  100 

c 

call  gsub(y,  q,  d,tn  +  t start) 
gcount  =  gcount  +  1 
go  to  100 
c 
c 
c 

entry  chemct  (tmk) 
c  - - 


c 

cd*  * 

3|C  %  4c 

* 

4c  4^  4c  ♦  ♦  ♦ 

4c4c4c*:tc4c4c4c4c*)fc******** 

cd 

cd 

chemct  (tmk) 

cd 

write  out  the 

values  of  the 

various  indicative  counters  that  the 

cd 

program  keeps , 

cd 

cd 

argument  list 

definition: 

cd 

tmk 

real 

a  floating  point  number  printed  i 

cd 

to  identify  the  call. 

cd 

cd 

output  variable 

definition: 

cd 

tmk 

real 

floating  point  identifier. 

cd 

gcount 

integer 

number  of  derivative  subroutine  calls 

cd 

since  the  last  call. 

cd 

rcount 

integer 

number  of  times  stepsize  was  reduced 

cd 

since  last  call. 

cd 

tgcnt 

integer 

total  of  gcount  to  this  call. 

cd 

trcnt 

integer 

total  of  rcount  to  this  call. 

cd 

cd*  ♦ 

4c  4c  4c  4c  4c  t  ♦ 

♦ 

4c  4c  *  4c  4c  * 

4c4c4c4c4c4c4c4c4c4c4c4c4;4c4c4c4c4c4c 

c 


r)3 


c 

c 


1000 


c 

c 


c 

c 

c 

c 

c 

cd*  ♦ 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 

cd 


tgcnt  -  tgcnt  +  gcount 
trcnt  =  trcnt  +  rcount 

print  out  indicative  counters. 

writedo,  1000)  tmk,  gcount,  rcount,  tgcnt, 
trcnt 

format (’  chemeq  indices;  tmk  =  ’,  lpelO.3, 

’  gcount,  rcount  =  ',  2i7,  »  totals;  ',  2i7) 

reset  counters . 
gcount  =  0 
rcount  =  0 
return 


entry  chemsp(epsmn,  epsmx,  dtmn,  tnot,  itermx,  ns,  ymn,  prt) 


chemspCepsmn,  epsmx,  dtmn,  tnot,  itermx,  prt) 


reset  any  local  control  parameters  if  their  respective  input 
values  are  greater  than  zero,  default  values  are  used  if  the 
input  values  are  zero  or  less  repectively. 


axgument 

list  definition 

epsmn 

real 

epsmx 

real 

dtmn 

real 

tnot 

real 

itermx 

i 

ns 

integer 

ymn(nd) 

real 

prt 

real 

the  maximum  relative  error  allowed  i 
for  convergence  of  the  corrector  step, 
default  value:  l.Oe-02 
this  number  provides  the  basis  for  i 
deciding  weather  convergence  can  be 
achieved  with  out  added  steps ize 
reduction,  if  eps/epsmin  is  greater 
than  epsmx  further  reduction  is 
applied, 

default  value  :  10.0 

the  smallest  steps ize  allowed.  i 

default  value:  l.Oe-15 

the  initial  value  of  the  independent  i 

variable  t, 

default  value:  0.0 

number  of  times  the  corrector  is  applied 
default  value:  1 

number  of  entries  in  ymin  to  reset  i 

minimum  values  allowed  for  y  i 

default  value:  l.Oe-20 
controls  the  output  of  chemsp.  any  i 
non  zero  value  suppresses  all  print 
output  from  this  entry. 
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cd 

c 

epsmin  =  l.Oe-02 

if(epsmn  .gt.  0.0) epsmin  =  epsmn 
if(epsmn  .gt.  0.0)sqreps  =  5 .0*sqrt (epsmin) 
epscl  =  1.0/epsmin 
epsmax  =  10.0 

if(epsrax  .gt.  0.0) epsmax  =  epsmx 

dtmin  =  l.Oe-15 

if (dtmn  .gt.  0.0) dtmin  =  dtmn 

tstart  =  tnot 

itermax  =  1 

if (itermx.gt .  0)  itermax  =  itermx 
do  i=l,ns 

ymin(i)  =  l.e-20 

if (ymn(i) .gt .0. )  ymin(i)  =  ymn(i) 
end  do 
c 

c  print  new  values  of  control  parameters. 
if(prt  .eq.  0.0)  then 

write(lo,  1001)  epsmn,  epsmx,  dtmn,  tnot,  itermx 
writedo,  1005)  ns 

if  (ns.gt.O)  writedo,  1006)  (ymin(i),  i=l,ns) 
end  if 

1001  format(’  initalize  **chemeq2'’  via  "chemsp"’,  /, 

’  epsmn,  epsmx,  dtmn,  tnot,  itermx  =  \  IpSglO.S) 

1005  format (’  ns  =  ’,15) 

1006  format (’  ymin:  ’,50el2,3) 

return 

c 

c 

c 

entry  chempr  (y,  n) 

c  - - 

c 

cd*  ********************************** 
cd 

cd  chempr  (y,  n) 

cd 

cd  chempr  may  be  called  whenever  an  error  occurs  that  can  be 

cd  attributed  to  the  results  of  chemeq.  a  partial  set  of  the  internal 

cd  variables  is  printed  as  a  diagnostic, 

cd 

cd  argument  list  definition: 

cd  y(n)  r  current  values  of  the  dependent  variable.  i 

cd  n  i  the  number  of  entries  in  y  and  ymin.  i 

cd 

cd*  ♦*♦**+**♦♦♦*************♦********* 


c 


write do,  1003)  dt,  tn,  dtmin 
do  45  i  =  l,n 

dtc  =  epsmin*y(i)/(abs(q(i)  -  d(i))  +  l.Oe-30) 

45  writeClo,  1004)  q(i).  d(i),  y(i),  rtau(i), 

^  q(i)-d(i),  ys(i).  y0(i),  yniin(i) 

return 

end 

subroutine  chemer 

c  - - 

c 

C  diagnostic  routine  for  stiff  o.d.e.  solver  -chemeq- 
print  1001 

1001  format(5(/),  »  library  version  of  -chemer-  called.’,  /, 

users  may  supply  their  own  version  for  diagnostics.’*,  /, 

’  no  arguments  are  required.’,  /, 

’program  will  continue  resetting  the  step  size  to  min-’,  /, 
.  ’imums  if  a  normal  return  is  made.’,  //, 

(stop  69)  executed  from  library  version  of  -chemer-’) 
c  ' 

stop  69 

end 
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10.2  Example  Driver  Code  and  Source  Term  Subroutines 


PROGRAM  CESIUM 

C  - 

C 

C  This  is  the  driver  program  for  the  seven-species  cesium 
C  mechanism  test  problem.  The  code  integrates  the  system 
C  MXCASE  times  using  differnt  values  of  the  chemeq2  variable 
C  epsmin  (set  by  passing  an  entry  from  array  EPS  through 
C  CHEMSP  before  each  integration) . 

C 

C  PROGRAM  SPECIFICATIONS. 

Q  - 

REAL  DSEC 

C 

REAL  Y(10),  YF(IO),  YMIN(IO),  YI(IO),  EPSIL(IO),  EPS(15) 

C 

INTEGER  SPSYM(7) 

C 

C  For  this  example,  the  external  subroutine  that  calculates  the 
C  source  terms  is  called  CSDFE, 

C 

EXTERNAL  CSDFE 
C 

DATA  YMIN/lO+l.OE-20/,  MXCASE/9/,  LO/16/ 

DATA  SPSYM/’02-’,  ’CS+>,  »CS’,  ’CS02»,  ’02’,  ’N2’,  ’NE’/ 

DATA  EPS/ 

,1,  .05,  .01,  ,005, 

.001,  .0005,  .0001,  .00005, 

.00001,  .000005,  .000001, 

5.e-7,  1 .e-7,5, e-8, 1 .e-8 
C  .  ,5.e-9,  l,e-9,5.e-10,l.e-10 

/ 

C 

1000  FORMAT(’CASE  NO,  ’,  15,  ’  PARAMETERS;’,  /, 

’  CONVERGENCE  PARAMETER  EPS  =  ’,  1PE10.3,  /, 

’  INNER  LOOP  LENGTH;’,  15) 

1001  F0RMAT(/,  ’  SPECIE  Y  -  INITAL  Y  -  FINAL  ’, 

’  Y  -  SOLUTION  REL  ERR’) 

1002  F0RMAT(5X,  A4,  1P3E1S.6,  E10.3) 

1003  FORMATC/,  ’  T  -  INITIAL  =  (’,  1PE10.3,  ’)  T  -  FINAL  =  (’, 

E10.3,  ’)’) 

1004  FORMATC/’  INTEGRATION  STATISTICS ;’ ) 

1005  FORMATC’  CPU  TIME  USED  FOR  INTEGRATION;’,  1PE10,3, 

’  SEC.,  CPU  TIME  NORMALIZED; ’ ,  18) 

1006  FORMATC’  SUM  OF  THE  RELATIVE  ERRORS  SQUARED;  ’,  1PE10.3) 

1007  FORMATC/) 

C 

C  Note  that  the  timing  routines  included  may  not  work  on 

C  all  systems.  Extra  timing  options  are  included  as  comments, 

C 

REAL+4  dtime,  delta,  tarray(2) 


integer  tnorm 
EXTERNAL  dtime 
delta  =  1. 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 


c 

c 


c 

c 


c 


INITIALIZE  CONTROL  PARAMETERS. 

TSCALE  IS  simply  a  normalization  factor  for  the  timing 
results.  It  can  be  used  to  compare  results  from  differnt 
machines  (by  setting  it  to  the  time  required  for  that 
machine  to  solve  a  standard  problem  of  some  sort)  or  to 
simply  make  the  timing  results  more  "friendly.” 

TSCALE  =  1.0/1024. 


INLP  allows  the  user  to  subdivide  the  interval  over  which 

each  test  is  run.  For  INLP=1,  CHEMEQZ  is  sent  the  full 

interval  TF-TI  (specified  below)  as  the  global  timestep. 
INLP  =1  o  r 


For  this  particular  test,  the  electron  number  density  is  not 
integrated.  The  other  five  reacting  species  are  integrated, 
^d  the  electron  density  is  found  through  charge  conservation. 
This  calculation  is  done  within  CSDFE.  Therefore,  NA  =  5  is 
the  number  of  equations  that  are  integrated,  but  NS  =  7  is  the 
number  of  species.  Species  to  be  integrated  must  be  placed  in 
first  NA  positions  within  the  Y  array.  CHEMEQ2  only  works  with 
these  first  NA  entries  since  NA  is  passed  in  the  argument  list 
e  ow,  ut  all  NS  values  are  available  to  and  used  by  CSDFE. 

NS  =  7 
NA  =  5 

"TI"  -  INITIAL  TIME,  "TF"  -  FINAL  TIME 
TI  =  0.0 
TF  =  1000.0 

DELTAT  =  (TF  -  TI)/INLP 

STORE  INITIAL(TI  =  0.0)  AND  FINALT(F  =  1000.0)  VALUES. 

02- 

YI(1)  =  6.200E+02 
YF(1)  =  2. 59 13949206 lD+04 

CS+ 

YI(2)  =  6.200E+02 
YF(2)  =  7.55718460300D+04 

CS 

YI(3)  =  l.OOOE+12 
YF(3)  =  1.53 19405 1722D+03 
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c 


CS02 


YI(4)  =  0. 

YF(4)  =  9.99999923516D+11 
C 

C  02 

YI(5)  =  3.600E+14 
YF(5)  =  3,590000000510+14 
C 

C  N2 

YI(6)  =  1.400E+15 
YF(6)  =  1.40000000000D+15 
C 

C  NE 

YI(7)  =  l.OOOE+02 
YF(7)  =  4,965789682390+04 


C 

C  LOOP  OVER  THE  TEST  CASES, 

00  30  ICASE  =  1 ,  MXCASE 

WRITE(LO,  1000)  ICASE,  EPS(ICASE),  INLP 
CALL  CHEMSP (EPS (ICASE),  0.,  0.,  TI,  5,  ns,  ymin,  0,) 
CPUT  =  0.0 
C 

C  RESET  "Y*‘  TO  INITIAL  VALUES  “YI". 

00  I  =  1,NS 
Y(I)  =  YI(I) 
end  do 
C 

C  SET  TIMER, 

C  Tl  =  SECN0S(0,0) 

delta  =  dtiine(t  array) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


c 

c 


INNER  LOOP  TO  OETERMINE  OVERHEAO  OR  RELATIVE  STARTING  EFFECIENCY 
OF  ITEGRATION  SCHEME  BEING  TESTEO, 

00  ISTEP  =  1,INLP 

CALL  INTEGRATOR. 

CALL  CHEMEq2(0ELTAT,  CSOFE,  NA,  Y) 


end  do 


CALCULATE  CPU  TIME  USEO  IN  THE  INTEGRATION  PROCESS, 
delta  -  dtime(t array) 

OSEC  =  SECNDS(Tl) 

OSEC  =  t array (1) 

OSEC  =  delta 
CPUT  =  CPUT  +  OSEC 
TNORM  =  INT(CPUT/TSCALE  +  .5) 

Calculate  final  electron  density  from  densities  of  other  charges  species 
Y(7)  =  Y(2)  -  Y(l) 
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C  CALCULATE  RELATIVE  ERROR. 

DO  I  =  1,NS 

EPSIL(I)  =  ABS(Y(I)  -  YF(I))/MIN(Y(I)  ,  YF(I)) 
end  do 
SUM  =  0.0 
DO  I  =  1,NS 

SUM  =  SUM  +  EPSIL(I)*+2 
end  do 

c  Root-mean-square  error  is  calculated  using  ns-1  (rather  than  ns) 
C  since  N2  is  inert. 


sum  -  sqrt(sum/real(ns-l)) 


PRINT  RESULTS. 

WRITE(L0,  1003)  TI,  TF 
WRITE(L0,  1001) 

DO  15  I  =  l.NS 

WRITE(L0,  1002)  SPSYM(I),  YI(I),  YF(I),  Y(I),  EPSIL(I) 
WRITE(L0,  1004) 

WRITE(L0,  1006)  SUM 
WRITE(L0,  1005)  CPUT,  TNQRM 
WRITE(*.699)  EPS(ICASE), 

CPUT,  TMORH,  sum 

format(lx,25HEPS,  time,  ticks,  error:  ,E7.1,2x  elO  4  2x 
fe  I5,2x,el0.4)  ’  ■  ’  ’ 

WRITE(L0,  1007) 

CALL  CHEMCT(TF) 

CONTINUE 
STOP  69 


subroutine  csdfe(y,  q,  d,  t) 

csdfe(y,  q,  d,  t) 
description: 

derivative  function  evaluator(gsub)  for  an  atmospheric  chemical 
relaxation  test  problem  involving  cesium  and  cesium  ions,  format¬ 
ion  and  loss  rates  are  calculated  for  this  set  of  "stiff  ordinary 
differential  equations"  that  was  suggested  by  by  d.  edelson  of 
bell  laboratories. 


argument  list  definitions: 


current  values  of  the  functions  plus  the  i/o 
extra  data  at  the  end  of  the  array  that  may  be 
passed  back  and  forth  between  "csdfe*'  and  the 
main  program,  locations  in  y(i)  which  represent 
the  functions  being  advanced  should  not  be 
tampered  with  here. 
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cd 

q(i) 

r 

total  formation  rates. 

i 

cd 

d(i) 

r 

total  loss  rates . 

i 

cd 

cd 

t 

r 

the  value  of  the  independent  variable. 

i 

cd  *********************************** 
c 

c  local  specifications. 


c  - 

real  ne ,  n2 

real  y(l),  q(l),  <i(i) 

c 

c  utilize  local  storage  for  varibles. 
o2m  =  y(l) 
csp  =  y(2) 
cs  =  y(3) 
cso2  =  y(4) 
o2  =  y(5) 
n2  =  y(6) 

c  write(63,*)  t 

c 

c  calculate  electron  density  for  local  use  and  transmission  back  to 
c  the  main  progrsim  via  y(7).  however  in  this  case  this  value  should 

c  not  be  trusted  since  "chemeq”  will  not  call  the  "gsub”  with  the 

c  latest  function  values  after  the  final  step  has  converged.  y(7) 

c  will  be  one  iteration  behind  in  this  case.  y(7)  and  y(6)  are 

c  examples  tho,  of  how  data  may  be  transfered  between  the  ”gsub*'  and 
c  the  main  program. 

ne  =  max(csp  -  o2m,  0.0) 
y(7)  =  ne 
c 

c  calculate  reaction  rates. 

crl  =  5-00e-08*o2m*csp 
cr2  =  1 .00e-12*csp*ne 
cr3  =  3.24e-03*cs 
cr4  =  4,00e-01*o2m 

cr5  =  1 .00e-31*o2*cs*(cs  +  cso2  +  n2  +  o2) 
cr6  =  1 .24e-30*o2*o2*ne 
cr7  1 . 00e-31*o2*n2*ne 

c  if (t .ge .700 , )  then 

c  cr4=  0 . 

c  cr6  =  0. 

c  cr7  =  0. 

c  end  if 

c 

c  calculate  total  formation  rates  (c(i))  and  total  loss  rates  (d(i)) 

c  for  each  species, 

c 

c  o2ro 
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q(l)  =  cr6  +  cr7 
<i(l)  =  crl  +  cr4 
c 

C  CS+ 

q(2)  =  cr3 

<i(2)  =  crl  +  cr2 
c 

c  cs 

q(3)  =  crl  +  cr2 
<i(3)  =  cr3  +  cr5 
c 

c  cso2 

q(4)  =  cr5 

C  q(4)  =  q(4)  -  1 .00e-31*o2*cs*cso2 

c  d(4)  =  -  1.00e-31*o2*cs*cso2 

c 

c  o2 

q(5)  =  crl  +  cr4 

d(5)  =  cr5  +  cr6  +  cr7 
c 

return 

end 
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10.3  Output  from  the  Sample  Programs 

Running  the  cesium  test  problem  as  given  with  an  compiler  flag  (or  equivalent)  will  result  in  the 

following  screen  output  describing  the  case,  unsealed  and  scaled  run  time,  and  the  resulting  rms  error: 


EPS,  time. 

ticks , 

error:  O.IE+OO 

0.9760E-03 

1 

0.4336E-02 

EPS,  time. 

ticks , 

error:  0.5E-01 

0.1952E-02 

2 

0.1035E-02 

EPS,  time. 

ticks. 

error:  0.1E~01 

0.5856E-02 

6 

0.1093E-03 

EPS,  time. 

ticks , 

error:  0.5E-02 

0.7808E-02 

8 

0.5876E-04 

EPS,  time. 

ticks , 

error:  0.1E‘“02 

0.2342E-01 

24 

0.9786E-05 

EPS,  time. 

ticks , 

error:  O.SE-OS 

0.3611E-01 

37 

0.4845E-05 

EPS,  time. 

ticks. 

error:  O.lE-03 

0.8491E-01 

87 

0.9995E“06 

EPS,  time. 

ticks. 

error:  0.5E-04 

0.1200E+00 

123 

0.5195E-06 

EPS,  time. 

ticks. 

error:  O.lE-04 

0.2538E+00 

260 

0.1154E-06 
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Of  course,  run  times  will  differ  on  different  platforms,  and  the  timing  routines  called  by  the  driver  routine 
may  not  be  available  on  all  systems. 

Additional  output  found  in  fort .  16  is  given  below.  This  file  holds  a  more  detailed  account  of  the  results 
of  each  integration,  including  a  count  of  the  number  of  times  the  source  term  subroutine  was  called  and  the 
number  of  timesteps  that  were  redone.  Included  below  is  this  information  for  the  last  calculation  in  the  test 
problem,  for  EPS  =  1.000E*“05. 

CASE  NO.  9  PARAMETERS; 

CONVERGENCE  PARAMETER  EPS  =  l.OOOE-05 

INNER  LOOP  LENGTH;  1 

initalize  *‘chemeq2’’  via  ‘^chemsp’’ 

epsmn,  epsmx,  dtmn,  tnot,  itermx  =  l.OOOE-OS  O.OOOE+00  O.OOOE+00  O.OOOE+00 
5 

ns  =  7 

ymin:  O.lOOE-19  O.lOOE-19  0.100E-i9  O.lOOE-19  O.lOOE-19 

O.lOOE-19  O.lOOE-19 

T  -  INITIAL  =  (  O.OOOE+00)  T  -  FINAL  =  (  l,000E+03) 


SPECIE 

Y  “  INITAL 

Y  -  FINAL 

Y  -  SOLUTION 

REL  ERR 

02“ 

5.200000E+02 

2.591395E+04 

2.591395E+04 

9 . 164E“08 

CS+ 

6.200000E+02 

7,557185E+04 

7.557185E+04 

9,755E“08 

CS 

l.OOOOOOE+12 

l,531941E+03 

1.531941E+03 

2.271E-07 

CS02 

O.OOOOOOE+00 

9.999999E+11 

9.999999E+11 

1.467E-08 

02 

3.600000E+14 

3.590000E+14 

3,590000E+14 

1.485E-09 

N2 

1,400000E+15 

1.400000E+15 

1.400000E+15 

O.OOOE+00 

NE 

l,000000E+02 

4.965790E+04 

4,965790E+04 

1.006E-07 

INTEGRATION  STATISTICS; 

SUM  OF  THE  RELATIVE  ERRORS  SQUARED;  l,154E-07 

CPU  TIME  USED  FOR  INTEGRATION;  2.538E-01  SEC.,  CPU  TIME  NORMALIZED; 
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chemeq  indices;  tmk  =  l.OOOE+03  gcount ,  rcount  =  149451  3  totals; 

313597  29 
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